{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "<div style=\"text-align: right\" align=\"right\"><i>Peter Norvig<br> Feb 2020<br>revised Jul 2025</i></div> \n",
    "\n",
    "# Solving Probability Problems with Monte Carlo Simulation\n",
    "\n",
    "There are at least three distinct approaches to discrete probability problems:\n",
    "- **Formula**: Apply a mathematical formula (perhaps involving [binomials](https://en.wikipedia.org/wiki/Binomial_distribution) or [combinations](https://en.wikipedia.org/wiki/Combination) or [Bayes' Rule](https://en.wikipedia.org/wiki/Bayes%27_theorem)).\n",
    "- **Enumeration**: Exhaustively list all possible outcomes and add up their probabilities. I cover this in [**another notebook**](http://nbviewer.jupyter.org/url/norvig.com/ipython/Probability.ipynb).\n",
    "- **Simulation**: Run a [Monte Carlo simulation](https://en.wikipedia.org/wiki/Monte_Carlo_method), in which, instead of considering all possible values for each random variable, we randomly select one outcome at each choice point, each one contingent on the previous choices, until we get a final outcome. Then we repeat the whole process *N* times to get a distribution of possible outcomes. \n",
    "\n",
    "# Pros of Monte Carlo Simulations\n",
    "\n",
    "Sometimes there is no simple mathematical formula, and it is impractical, inconvenient, or even impossible to explicitly enumerate all possible outcomes. Then a simulation may be just the thing, because it is:\n",
    "\n",
    "- **Practical**: Suppose some problem has a state space with, say, 27 variables, each with 1,000 possible values. We can't enumerate all the states on our computer, because there are more states than atoms in the universe. But we could do a simulation where we make a choice for each of the 27 variables, and repeat a million times. That's practical, and can give a good estimate of the true distribution.\n",
    "\n",
    "- **Convenient**: The code for a simulation is simpler. You can keep track of the state of the world in ordinary programming variables, and update them with `if` and assignment statements. The code for an enumeration has to package up the state of the world into a structure, and keep track of every possible update of every possible state. Thus, a simple assignment statement in the simulation because a complex nested loop in the enumeration.\n",
    "\n",
    "- **Possible**: Consider a game like rock-paper-scissors: it is impossible to fully enumerate the possible games, because a game can last forever (when both players continually make the same play). In a simulation an infinite game is theoretically possible, but in practice all simulated games will end fairly quickly. Or consider a game where there is an infinite number of possible moves on each turn (i.e., you can pick a real number, not just an integer). Yes, we could do tricks to make thee games finite. But in a simulation we don't have to worry about it.\n",
    " \n",
    "# Cons of Monte Carlo Simulations\n",
    "\n",
    "- **Inexactness**: A simulation gives only approximate answers. The chance of a fair coin flip coming up heads is (by definition) exactly 50%, but a simulation will in general give a slightly different answer each run, with an error that is inversely proportional to the square root of the number of trials (meaning that you need to do 100 times more trials to get 10 times more accuracy). \n",
    "\n",
    "- **Imbalanced data**: A problem arises when there are situations that are very rare, but very consequential.  Suppose we are simulating the flight of an airplane or the running of a power plant, and there is a very small possibility that a sequence of events could lead to a catastrophic failure. We might have to run billions of trials to see this even once, and we'd like to see it many times to get a good estimate of its probability.\n",
    "\n",
    "\n",
    "\n",
    "# Simulating a Game: Monopoly\n",
    "\n",
    "<img src=\"http://buckwolf.org/a.abcnews.com/images/Entertainment/ho_hop_go_050111_t.jpg\" style=\"float:right\">\n",
    "\n",
    "Consider [problem 84](https://projecteuler.net/problem=84) from the excellent [Project Euler](https://projecteuler.net), which asks: in the board game Monopoly, what is the probability of ending a roll of the dice on each square on the board? Not every square will have the same probability, because chance and community chest cards direct players to specific squares, and there are multiple ways to end up on the \"Jail\" square.  This looks like a problem where a simulation would be a good approach. I'll define the function `monopoly(N)` to simulate a game for `N` random dice rolls and generate the squares where the player ends each roll.\n",
    "\n",
    "I will represent the board by a list of 40 names of properties, using Project Euler's abbreviations. Squares  are represented by indexes into the board–integers from 0 to 39.    I'll implement a deck of cards (chance or community chest) as an infinite cycle of shuffled cards. Each card is either a square number (meaning go to that square), the number -3 (meaning go back three spaces), a set of square numbers (either the set of utilities or the set of railroads, meaning that we advance to the nearest utility/railroad), or `None`, meaning stay at the same square (but gain or lose some money, which this simulation does not track).\n",
    "\n",
    "\n",
    "Here's the code for the simulation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import random\n",
    "from collections import Counter\n",
    "from itertools   import combinations, cycle, count as count_from\n",
    "from statistics  import mean\n",
    "from typing      import List, Iterable, Iterator, Callable, Tuple"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The Monopoly board, as specified by https://projecteuler.net/problem=84\n",
    "\n",
    "Square  = int       # Data type for a square on the board: the square number\n",
    "Board   = List[str] # Data type for board; a list of square names (not numbers)\n",
    "Deck    = Iterator  # Data type for a deck of cards; you can take the `next` one.\n",
    "\n",
    "board: Board = \"\"\"\n",
    " GO    A1  CC1  A2   T1  R1  B1   CH1  B2  B3 \n",
    " JAIL  C1  U1   C2   C3  R2  D1   CC2  D2  D3  \n",
    " FP    E1  CH2  E2   E3  R3  F1   F2   U2  F3  \n",
    " G2J   G1  G2   CC3  G3  R4  CH3  H1   T2  H2\"\"\".split()    \n",
    "\n",
    "globals().update({board[s]: s for s in range(len(board))}) # Square names are global variables\n",
    "\n",
    "RRs = {R1, R2, R3, R4} # The 4 railroads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Chance and Community Chest card decks\n",
    "\n",
    "def deck(cards: list) -> Deck:\n",
    "    \"\"\"Make a deck of cards: an infinite iterable cycling through the (shuffled) cards.\"\"\"\n",
    "    random.shuffle(cards)\n",
    "    return cycle(cards)\n",
    "    \n",
    "CC_cards = deck([None] * 14 + [GO, JAIL])\n",
    "CH_cards = deck([None] *  6 + [GO, JAIL, C1, E3, H2, R1, RRs, RRs, {U1, U2}, -3])\n",
    "\n",
    "def do_card(card, square: Square) -> Square:\n",
    "    \"\"\"Update location from `square` to new location, depending on what `card` says.\"\"\"\n",
    "    return (     square     if card is None                    # Don't move (card is about money)\n",
    "            else square - 3 if card == -3                      # Go back 3 spaces\n",
    "            else card       if isinstance(card, Square)        # Go to square named on card\n",
    "            else min([s for s in card if s > square] or card)) # Advance to nearest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The simulation\n",
    "\n",
    "def monopoly(N: int) -> Iterable[Square]:\n",
    "    \"\"\"Simulate a Monopoly game for `N` dice rolls, yielding the squares visited.\"\"\"\n",
    "    square = GO                                      # Current location on board\n",
    "    doubles = 0                                      # Count of consecutive doubles rolled\n",
    "    for _ in range(N):\n",
    "        d1, d2 = roll(), roll()                      # Roll two dice\n",
    "        square = (square + d1 + d2) % len(board)     # Move ahead, maybe pass Go\n",
    "        doubles = (doubles + 1 if d1 == d2 else 0)\n",
    "        if square == G2J or doubles == 3: \n",
    "            square = JAIL                            # Go to Jail\n",
    "            doubles = 0\n",
    "        elif board[square].startswith('CC'):         # Community Chest card\n",
    "            square = do_card(next(CC_cards), square)\n",
    "        elif board[square].startswith('CH'):         # Chance card\n",
    "            square = do_card(next(CH_cards), square)\n",
    "        yield square\n",
    "\n",
    "def roll(d=6) -> int: \n",
    "    \"\"\"Random roll of a d-sided die.\"\"\"\n",
    "    return random.randint(1, d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's run the simulation for a million rolls and look at the resulting counts of square visits:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Counter({9: 23042,\n",
       "         17: 25898,\n",
       "         22: 10376,\n",
       "         27: 26773,\n",
       "         35: 24422,\n",
       "         3: 21704,\n",
       "         19: 30994,\n",
       "         26: 27071,\n",
       "         32: 26113,\n",
       "         12: 26086,\n",
       "         10: 62259,\n",
       "         15: 29294,\n",
       "         24: 31746,\n",
       "         31: 26721,\n",
       "         0: 31009,\n",
       "         4: 23348,\n",
       "         13: 23671,\n",
       "         21: 28515,\n",
       "         14: 24621,\n",
       "         23: 27259,\n",
       "         5: 29258,\n",
       "         36: 8789,\n",
       "         6: 22708,\n",
       "         11: 26968,\n",
       "         38: 21831,\n",
       "         28: 27853,\n",
       "         33: 23602,\n",
       "         39: 26357,\n",
       "         1: 21215,\n",
       "         18: 29451,\n",
       "         20: 28791,\n",
       "         16: 27738,\n",
       "         29: 26122,\n",
       "         25: 30728,\n",
       "         34: 25102,\n",
       "         2: 19087,\n",
       "         8: 23233,\n",
       "         37: 21658,\n",
       "         7: 8587})"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P = Counter(monopoly(1_000_000))\n",
    "P"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That will be easier to visualize as a bar chart, with square names on the x-axis. On the y-axis I want probabilities, not counts, so I will introduce the class `ProbDist` for \"Probability Distribution.\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "Probability = float # Data Type for a probability, a float between 0.0 and 1.0, inclusive.\n",
    "\n",
    "class ProbDist(Counter):\n",
    "    \"\"\"Probability Distribution. Like a Counter, but normalized so that values sum to 1.0.\"\"\"\n",
    "    def __init__(self, outcomes=()):\n",
    "        super().__init__(outcomes)\n",
    "        N = sum(self.values())\n",
    "        for x in self:\n",
    "            self[x] /= N\n",
    "\n",
    "def monopoly_barchart(counts: Counter, board=board):\n",
    "    \"\"\"A bar chart of {square: probability} with names for squares.\"\"\"\n",
    "    plt.figure(figsize=[11, 7])\n",
    "    plt.xticks(rotation=90)\n",
    "    N = len(board)\n",
    "    plt.plot([0, N - 1], [1/N, 1/N], 'r:') # Average value\n",
    "    P_barchart(ProbDist({board[i]: counts[i] for i in range(len(board))}))\n",
    "\n",
    "def P_barchart(P: ProbDist): \n",
    "    \"\"\"A bar chart of a discrete probability distribution.\"\"\" \n",
    "    plt.xlabel(\"Outcomes\")   \n",
    "    plt.ylabel(\"Probability\")\n",
    "    plt.grid(axis='y', ls=':')\n",
    "    plt.bar(P.keys(), P.values())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6QAAAJpCAYAAAC6vv8DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB9GUlEQVR4nO3deXxTVf7/8fdNuoFQdoowZd9lHBDUqQqICsgyLgMjOo4Lgg7iCqOOLOqIIu4gyvIdFXcdf444g1BZFAUUHARBXMoilEVpVbYWQbok5/dHTEjapE3SW24pr+fjkUdP7705PZ+eJfeT3CSWMcYIAAAAAIBjzOV0AwAAAAAAJyYSUgAAAACAI0hIAQAAAACOICEFAAAAADiChBQAAAAA4AgSUgAAAACAI0hIAQAAAACOICEFAAAAADgiwekGVEVer1e7d+9W7dq1ZVmW080BAAAAgOOKMUYHDx5U06ZN5XJFfh2UhDSM3bt3Kz093elmAAAAAMBxbdeuXfrNb34TcT8JaRi1a9eW5PvnpaamOtwaAAAAADi+5OfnKz09PZBbRUJCGob/Mt3U1FQSUgAAAACIU3lvgeRDjQAAAAAAjiAhBQAAAAA4goQUAAAAAOAIElIAAAAAgCNISAEAAAAAjiAhBQAAAAA4goQUAAAAAOAIElIAAAAAgCNISAEAAAAAjiAhBQAAAAA4goQUAAAAAOAIElIAAAAAgCNISAEAAAAAjiAhBQAAAAA4goQUAAAAAOAIElIAAAAAgCNISAEAAAAAjiAhBQAAAAA4goQUAAAAAOAIElIAAAAAgCNISAEAAAAAjiAhBQAAAAA4IsHpBgBAvFrevSDu+25/eJCNLQEAAEA8eIUUAAAAAOAIElIAAAAAgCNISAEAAAAAjiAhBQAAAAA4goQUAAAAAOAIElIAAAAAgCNISAEAAAAAjiAhBQAAAAA4goQUAAAAAOAIElIAAAAAgCNISAEAAAAAjiAhBQAAAAA4goQUAAAAAOAIElIAAAAAgCMcT0hnzpypVq1aKSUlRd27d9eKFSvKPH7ZsmXq3r27UlJS1Lp1a82ePbvUMQcOHNBNN92kk08+WSkpKerUqZMyMzMrKwQAAAAAQBwcTUjffPNN3X777ZowYYLWrVunnj17asCAAdq5c2fY47OzszVw4ED17NlT69at0/jx43Xrrbfq7bffDhxTWFiovn37avv27fr3v/+tTZs26dlnn1WzZs2OVVgAAAAAgChYxhjj1B8/88wzddppp2nWrFmBbZ06ddIll1yiKVOmlDr+73//u+bNm6esrKzAtlGjRumLL77QqlWrJEmzZ8/WY489po0bNyoxMTGuduXn56tOnTrKy8tTampqXHUAqHwt714Q9323PzzIxpYAAAAgWLQ5lWOvkBYWFmrt2rXq169fyPZ+/fpp5cqVYe+zatWqUsf3799fa9asUVFRkSRp3rx5ysjI0E033aS0tDR16dJFDz30kDweT8S2FBQUKD8/P+QmKXAfj8cTtlxcXBxS9nq9ZZaLiopCyv7nAvxlY0ypsqSQstfrDSkXFxeXWfZ4PCFlYiKm6hSTS744kt1GVqmyUbLb99MKlBVSrooxVcd+IiZiIiZiIiZiIqYTN6byOJaQ7tmzRx6PR2lpaSHb09LSlJubG/Y+ubm5YY8vLi7Wnj17JEnbtm3Tv//9b3k8HmVmZmrixIl64oknNHny5IhtmTJliurUqRO4paenS5K++uorSVJWVlbgVdkNGzZoy5YtkqR169YpOztbkrR69Wrt2rVLkrRy5Url5ORIkpYvXx5o29KlS3XgwAFJ0uLFi3Xw4EFJUmZmpo4cOaLi4mJlZmaquLhYR44cCbzv9eDBg1q8eLEk3/tjly5dGvgfLl++XJKUk5MTSOR37dql1atXS/Jd5rxu3TpJ0pYtW7RhwwZiIqZqE1PzWpIkTTrNo7QavvKjZ3iUmiQlu33lZLeUmuQrS1JaDd/xVTWm6thPxERMxERMxERMxHTixbRp0yZFw7FLdnfv3q1mzZpp5cqVysjICGyfPHmyXnnlFW3cuLHUfdq3b6/hw4dr3LhxgW2ffPKJzjnnHOXk5KhJkyZq3769jhw5ouzsbLndbknSk08+qcceeyzQcSUVFBSooKAg8Ht+fr7S09O1b98+1atXL5Ddu93ukHJxcbEsywqUXS6XXC5XxHJRUZHcbnegnJCQIMuyAmXJ9yxHcDkxMVHGmEDZ6/XK4/EEyl6vVwkJCRHLHo9HxphAOVwcxERMx2tMHe5ZJK8sJbuNCj2+10KPln1JaYFHsiQluaUCjyVLRkluadPkwVUypurYT8RETMRETMRETMR04sWUn5+v+vXrl3vJrmMJaWFhoWrWrKm33npLl156aWD7bbfdpvXr12vZsmWl7tOrVy9169ZNTz31VGDbO++8o8suu0yHDx9WYmKievfurcTERL3//vuBY9577z0NHDhQBQUFSkpKKrdtvIcUOD7wHlIAAICqqcq/hzQpKUndu3fXkiVLQrYvWbJEZ511Vtj7ZGRklDp+8eLF6tGjR+ADjM4++2x9++23geunJWnz5s06+eSTo0pGAQAAAADHhqNf+zJ27Fg999xzmjNnjrKysjRmzBjt3LlTo0aNkiSNGzdOV199deD4UaNGaceOHRo7dqyysrI0Z84cPf/887rjjjsCx9x4443au3evbrvtNm3evFkLFizQQw89pJtuuumYxwcAAAAAiCzByT8+bNgw7d27V5MmTVJOTo66dOmizMxMtWjRQpLvzbPB30naqlUrZWZmasyYMZoxY4aaNm2q6dOna8iQIYFj0tPTtXjxYo0ZM0annnqqmjVrpttuu01///vfj3l8AAAAAIDIHP0e0qqK95ACxwfeQwoAAFA1Vfn3kAIAAAAATmwkpAAAAAAAR5CQAgAAAAAcQUIKAAAAAHAECSkAAAAAwBEkpAAAAAAAR5CQAgAAAAAcQUIKAAAAAHAECSkAAAAAwBEkpAAAAAAAR5CQAgAAAAAcQUIKAAAAAHAECSkAAAAAwBEkpAAAAAAAR5CQAgAAAAAcQUIKAAAAAHAECSkAAAAAwBEkpAAAAAAAR5CQAgAAAAAcQUIKAAAAAHAECSkAAAAAwBEkpAAAAAAAR5CQAgAAAAAcQUIKAAAAAHAECSkAAAAAwBEkpAAAAAAAR5CQAgAAAAAcQUIKAAAAAHAECSkAAAAAwBEkpAAAAAAAR5CQAgAAAAAcQUIKAAAAAHAECSkAAAAAwBEkpAAAAAAAR5CQAgAAAAAcQUIKAAAAAHAECSkAAAAAwBEkpAAAAAAAR5CQAgAAAAAcQUIKAAAAAHAECSkAAAAAwBEkpAAAAAAAR5CQAgAAAAAcQUIKAAAAAHAECSkAAAAAwBEkpAAAAAAAR5CQAgAAAAAcQUIKAAAAAHAECSkAAAAAwBEkpAAAAAAAR5CQAgAAAAAcQUIKAAAAAHAECSkAAAAAwBEkpAAAAAAAR5CQAgAAAAAcQUIKAAAAAHAECSkAAAAAwBEkpAAAAAAAR5CQAgAAAAAcQUIKAAAAAHAECSkAAAAAwBEkpAAAAAAAR5CQAgAAAAAcQUIKAAAAAHAECSkAAAAAwBEkpAAAAAAAR5CQAgAAAAAcQUIKAAAAAHAECSkAAAAAwBEkpAAAAAAAR5CQAgAAAAAc4XhCOnPmTLVq1UopKSnq3r27VqxYUebxy5YtU/fu3ZWSkqLWrVtr9uzZIftffPFFWZZV6nbkyJHKDAMAAAAAECNHE9I333xTt99+uyZMmKB169apZ8+eGjBggHbu3Bn2+OzsbA0cOFA9e/bUunXrNH78eN166616++23Q45LTU1VTk5OyC0lJeVYhAQAAAAAiFKCk3/8ySef1IgRIzRy5EhJ0rRp07Ro0SLNmjVLU6ZMKXX87Nmz1bx5c02bNk2S1KlTJ61Zs0aPP/64hgwZEjjOsiw1adLkmMQAAAAAAIiPY6+QFhYWau3aterXr1/I9n79+mnlypVh77Nq1apSx/fv319r1qxRUVFRYNvPP/+sFi1a6De/+Y0GDx6sdevWldmWgoIC5efnh9wkyePxBH6GKxcXF4eUvV5vmeWioqKQsjEmpGyMKVWWFFL2er0h5eLi4jLLHo8npExMxFSdYnLJF0ey28gqVTZKdvt+WoGyQspVMabq2E/EREzEREzEREzEdOLGVB7HEtI9e/bI4/EoLS0tZHtaWppyc3PD3ic3Nzfs8cXFxdqzZ48kqWPHjnrxxRc1b948vfHGG0pJSdHZZ5+tLVu2RGzLlClTVKdOncAtPT1dkvTVV19JkrKyspSVlSVJ2rBhQ6CudevWKTs7W5K0evVq7dq1S5K0cuVK5eTkSJKWL18eaNvSpUt14MABSdLixYt18OBBSVJmZqaOHDmi4uJiZWZmqri4WEeOHFFmZqYk6eDBg1q8eLEk6cCBA1q6dGngf7h8+XJJUk5OTiCR37Vrl1avXi3Jd5mzPyHfsmWLNmzYQEzEVG1ial5LkqRJp3mUVsNXfvQMj1KTpGS3r5zsllKTfGVJSqvhO76qxlQd+4mYiImYiImYiImYTryYNm3apGhYxp86H2O7d+9Ws2bNtHLlSmVkZAS2T548Wa+88oo2btxY6j7t27fX8OHDNW7cuMC2Tz75ROecc45ycnLCXqbr9Xp12mmnqVevXpo+fXrYthQUFKigoCDwe35+vtLT07Vv3z7Vq1cvkN273e6QcnFxsSzLCpRdLpdcLlfEclFRkdxud6CckJAgy7ICZcn3LEdwOTExUcaYQNnr9crj8QTKXq9XCQkJEcsej0fGmEA5XBzEREzHa0wd7lkkrywlu40KPb7XQo+WfUlpgUeyJCW5pQKPJUtGSW5p0+TBVTKm6thPxERMxERMxERMxHTixZSfn6/69esrLy9PqampisSxhLSwsFA1a9bUW2+9pUsvvTSw/bbbbtP69eu1bNmyUvfp1auXunXrpqeeeiqw7Z133tFll12mw4cPKzExMezfuv766/Xdd9/pvffei6pt+fn5qlOnTrn/PADOann3grjvu/3hQTa2BAAAAMGizakcu2Q3KSlJ3bt315IlS0K2L1myRGeddVbY+2RkZJQ6fvHixerRo0fEZNQYo/Xr1+vkk0+2p+EAAAAAAFs4+rUvY8eO1XPPPac5c+YoKytLY8aM0c6dOzVq1ChJ0rhx43T11VcHjh81apR27NihsWPHKisrS3PmzNHzzz+vO+64I3DM/fffr0WLFmnbtm1av369RowYofXr1wfqBAAAAABUDY5+7cuwYcO0d+9eTZo0STk5OerSpYsyMzPVokULSb43zwZ/J2mrVq2UmZmpMWPGaMaMGWratKmmT58e8pUvBw4c0A033KDc3FzVqVNH3bp10/Lly3XGGWcc8/gAAAAAAJE59h7Sqoz3kALHB95DCgAAUDVV+feQAgAAAABObCSkAAAAAABHkJACAAAAABxBQgoAAAAAcAQJKQAAAADAESSkAAAAAABHkJACAAAAABxBQgoAAAAAcAQJKQAAAADAESSkAAAAAABHkJACAAAAABxBQgoAAAAAcAQJKQAAAADAESSkAAAAAABHkJACAAAAABxBQgoAAAAAcAQJKQAAAADAESSkAAAAAABHkJACAAAAABxBQgoAAAAAcAQJKQAAAADAESSkAAAAAABHkJACAAAAABxBQgoAAAAAcAQJKQAAAADAESSkAAAAAABHkJACAAAAABxBQgoAAAAAcAQJKQAAAADAESSkAAAAAABHkJACAAAAABxBQgoAAAAAcAQJKQAAAADAESSkAAAAAABHkJACAAAAABxBQgoAAAAAcAQJKQAAAADAESSkAAAAAABHkJACAAAAABxBQgoAAAAAcAQJKQAAAADAESSkAAAAAABHkJACAAAAABxBQgoAAAAAcAQJKQAAAADAESSkAAAAAABHkJACAAAAABxBQgoAAAAAcAQJKQAAAADAESSkAAAAAABHkJACAAAAABxBQgoAAAAAcAQJKQAAAADAESSkAAAAAABHkJACAAAAABxBQgoAAAAAcAQJKQAAAADAESSkAAAAAABHkJACAAAAABxBQgoAAAAAcAQJKQAAAADAESSkAAAAAABHkJACAAAAABxBQgoAAAAAcAQJKQAAAADAESSkAAAAAABHkJACAAAAABxBQgoAAAAAcAQJKQAAAADAESSkAAAAAABHkJACAAAAABxBQgoAAAAAcITjCenMmTPVqlUrpaSkqHv37lqxYkWZxy9btkzdu3dXSkqKWrdurdmzZ0c89l//+pcsy9Ill1xic6sBAAAAABXlaEL65ptv6vbbb9eECRO0bt069ezZUwMGDNDOnTvDHp+dna2BAweqZ8+eWrduncaPH69bb71Vb7/9dqljd+zYoTvuuEM9e/as7DAAAAAAAHFwNCF98sknNWLECI0cOVKdOnXStGnTlJ6erlmzZoU9fvbs2WrevLmmTZumTp06aeTIkbruuuv0+OOPhxzn8Xh05ZVX6v7771fr1q2PRSgAAAAAgBg5lpAWFhZq7dq16tevX8j2fv36aeXKlWHvs2rVqlLH9+/fX2vWrFFRUVFg26RJk9SoUSONGDEiqrYUFBQoPz8/5Cb5Elv/z3Dl4uLikLLX6y2zXFRUFFI2xoSUjTGlypJCyl6vN6RcXFxcZtnj8YSUiYmYqlNMLvniSHYbWaXKRslu308rUFZIuSrGVB37iZiIiZiIiZiIiZhO3JjK41hCumfPHnk8HqWlpYVsT0tLU25ubtj75Obmhj2+uLhYe/bskSR98sknev755/Xss89G3ZYpU6aoTp06gVt6erok6auvvpIkZWVlKSsrS5K0YcMGbdmyRZK0bt06ZWdnS5JWr16tXbt2SZJWrlypnJwcSdLy5csDbVu6dKkOHDggSVq8eLEOHjwoScrMzNSRI0dUXFyszMxMFRcX68iRI8rMzJQkHTx4UIsXL5YkHThwQEuXLg38D5cvXy5JysnJCSTyu3bt0urVqyX5LnNet26dJGnLli3asGEDMRFTtYmpeS1Jkiad5lFaDV/50TM8Sk2Skt2+crJbSk3ylSUprYbv+KoaU3XsJ2IiJmIiJmIiJmI68WLatGmTomEZf+p8jO3evVvNmjXTypUrlZGREdg+efJkvfLKK9q4cWOp+7Rv317Dhw/XuHHjAts++eQTnXPOOcrJydFJJ52kU089VTNnztSAAQMkSddee60OHDig//znPxHbUlBQoIKCgsDv+fn5Sk9P1759+1SvXr1Adu92u0PKxcXFsiwrUHa5XHK5XBHLRUVFcrvdgXJCQoIsywqUJd+zHMHlxMREGWMCZa/XK4/HEyh7vV4lJCRELHs8HhljAuVwcRATMR2vMXW4Z5G8spTsNir0+F4LPVr2JaUFHsmSlOSWCjyWLBkluaVNkwdXyZiqYz8REzEREzEREzER04kXU35+vurXr6+8vDylpqYqEscS0sLCQtWsWVNvvfWWLr300sD22267TevXr9eyZctK3adXr17q1q2bnnrqqcC2d955R5dddpkOHz6sr7/+Wt26dZPb7Q7s979s7XK5tGnTJrVp06bctuXn56tOnTrl/vMAOKvl3Qvivu/2hwfZ2BIAAAAEizancuyS3aSkJHXv3l1LliwJ2b5kyRKdddZZYe+TkZFR6vjFixerR48eSkxMVMeOHfXll19q/fr1gdtFF12kPn36aP369YFLcQEAAAAAzktw8o+PHTtWV111lXr06KGMjAz985//1M6dOzVq1ChJ0rhx4/T999/r5ZdfliSNGjVKzzzzjMaOHavrr79eq1at0vPPP6833nhDkpSSkqIuXbqE/I26detKUqntAAAAAABnOZqQDhs2THv37tWkSZOUk5OjLl26KDMzUy1atJDke/Ns8HeStmrVSpmZmRozZoxmzJihpk2bavr06RoyZIhTIQAAAAAA4uTYe0irMt5DChwfeA8pAABA1VTl30MKAAAAADixkZACAAAAABxBQgoAAAAAcAQJKQAAAADAESSkAAAAAABHkJACAAAAABxBQgoAAAAAcAQJKQAAAADAESSkAAAAAABHkJACAAAAABxBQgoAAAAAcAQJKQAAAADAESSkAAAAAABHkJACAAAAABxBQgoAAAAAcAQJKQAAAADAESSkAAAAAABHkJACAAAAABxBQgoAAAAAcAQJKQAAAADAEXElpC+++KIOHz5sd1sAAAAAACeQuBLScePGqUmTJhoxYoRWrlxpd5sAAAAAACeAuBLS7777Tq+++qr279+vPn36qGPHjnrkkUeUm5trd/sAAAAAANVUXAmp2+3WRRddpLlz52rXrl264YYb9Nprr6l58+a66KKL9N///lder9futgIAAAAAqpEKf6hR48aNdfbZZysjI0Mul0tffvmlrr32WrVp00YfffSRDU0EAAAAAFRHcSekP/zwgx5//HGdcsopOvfcc5Wfn6/58+crOztbu3fv1h//+Eddc801drYVAAAAAFCNJMRzpz/84Q9atGiR2rdvr+uvv15XX3216tevH9hfo0YN/e1vf9PUqVNtaygAAAAAoHqJKyFt3Lixli1bpoyMjIjHnHzyycrOzo67YQAAAACA6i2uS3Z79+6t0047rdT2wsJCvfzyy5Iky7LUokWLirUOAAAAAFBtxZWQDh8+XHl5eaW2Hzx4UMOHD69wowAAAAAA1V9cCakxRpZlldr+3XffqU6dOhVuFAAAAACg+ovpPaTdunWTZVmyLEvnn3++EhKO3t3j8Sg7O1sXXnih7Y0EAAAAAFQ/MSWkl1xyiSRp/fr16t+/v2rVqhXYl5SUpJYtW2rIkCG2NhAAAAAAUD3FlJDed999kqSWLVtq2LBhSklJqZRGAQAAAACqv7i+9uWaa66xux0AAAAAgBNM1Alp/fr1tXnzZjVs2FD16tUL+6FGfvv27bOlcQAAAACA6ivqhHTq1KmqXbt2oFxWQgoAAAAAQHmiTkiDL9O99tprK6MtAAAAAIATSNQJaX5+ftSVpqamxtUYAAAAAMCJI+qEtG7duuVepmuMkWVZ8ng8FW4YAAAAAKB6izoh/fDDDyuzHQAAAACAE0zUCWnv3r0rsx0AAAAAgBNM1Anphg0b1KVLF7lcLm3YsKHMY0899dQKNwwAAAAAUL1FnZB27dpVubm5aty4sbp27SrLsmSMKXUc7yEFAAAAAEQj6oQ0OztbjRo1CpQBAAAAAKiIqBPSFi1ahC0DAAAAABCPqBPSkjZt2qSnn35aWVlZsixLHTt21C233KIOHTrY2T4AAAAAQDXliudO//73v9WlSxetXbtWv/vd73Tqqafq888/V5cuXfTWW2/Z3UYAAAAAQDUU1yukd911l8aNG6dJkyaFbL/vvvv097//XX/6059saRwAAAAAoPqK6xXS3NxcXX311aW2/+Uvf1Fubm6FGwUAAAAAqP7iSkjPPfdcrVixotT2jz/+WD179qxwowAAAAAA1V/Ul+zOmzcvUL7ooov097//XWvXrtXvf/97SdKnn36qt956S/fff7/9rQQAAAAAVDuWMcZEc6DLFd2LqZZlyePxVKhRTsvPz1edOnWUl5en1NRUp5sDIIKWdy+I+77bHx5kY0sAAAAQLNqcKupXSL1ery0NAwAAAABAivM9pAAAAAAAVFRcX/siSYcOHdKyZcu0c+dOFRYWhuy79dZbK9wwAAAAAED1FldCum7dOg0cOFCHDx/WoUOHVL9+fe3Zs0c1a9ZU48aNSUgBAAAAAOWK65LdMWPG6A9/+IP27dunGjVq6NNPP9WOHTvUvXt3Pf7443a3EQAAAABQDcWVkK5fv15/+9vf5Ha75Xa7VVBQoPT0dD366KMaP3683W0EAAAAAFRDcSWkiYmJsixLkpSWlqadO3dKkurUqRMoAwAAAABQlrjeQ9qtWzetWbNG7du3V58+fXTvvfdqz549euWVV/Tb3/7W7jYCAAAAAKqhuBLShx56SAcPHpQkPfDAA7rmmmt04403qm3btnrhhRdsbSAAAIDdWt69oEL33/7wIJtaAgAntrgS0h49egTKjRo1UmZmpm0NAgAAAACcGOL+HlJJ+vHHH7Vp0yZZlqUOHTqoUaNGdrULAAAAAFDNxfWhRvn5+brqqqvUrFkz9e7dW7169VLTpk31l7/8RXl5eXa3EQAAAABQDcWVkI4cOVL/+9//NH/+fB04cEB5eXmaP3++1qxZo+uvv97uNgIAAAAAqqG4LtldsGCBFi1apHPOOSewrX///nr22Wd14YUX2tY4AAAAAED1FVdC2qBBA9WpU6fU9jp16qhevXoVbhTKxicDAsDxjXUcAACfuC7ZnThxosaOHaucnJzAttzcXN1555265557bGscAAAAAKD6ivoV0m7dusmyrMDvW7ZsUYsWLdS8eXNJ0s6dO5WcnKyffvpJf/3rX+1vKQAAQBVVkVe9ecUbwIks6oT0kksuqcRmAAAAAABONFEnpPfdd1+lNGDmzJl67LHHlJOTo1NOOUXTpk1Tz549Ix6/bNkyjR07Vl9//bWaNm2qu+66S6NGjQrsnzt3rh566CF9++23KioqUrt27fS3v/1NV111VaW0HwAAAAAQn7g+1Mhv7dq1ysrKkmVZ6ty5s7p16xbT/d98803dfvvtmjlzps4++2z93//9nwYMGKBvvvkmcClwsOzsbA0cOFDXX3+9Xn31VX3yyScaPXq0GjVqpCFDhkiS6tevrwkTJqhjx45KSkrS/PnzNXz4cDVu3Fj9+/evSLgAAAAAABvFlZD++OOPuvzyy/XRRx+pbt26MsYoLy9Pffr00b/+9S81atQoqnqefPJJjRgxQiNHjpQkTZs2TYsWLdKsWbM0ZcqUUsfPnj1bzZs317Rp0yRJnTp10po1a/T4448HEtJzzz035D633XabXnrpJX388cckpAAAAABQhcT1Kbu33HKL8vPz9fXXX2vfvn3av3+/vvrqK+Xn5+vWW2+Nqo7CwkKtXbtW/fr1C9ner18/rVy5Mux9Vq1aVer4/v37a82aNSoqKip1vDFGH3zwgTZt2qRevXpFbEtBQYHy8/NDbpLk8XgCP8OVi4uLQ8per7fMclFRUUjZGBNSNsaUKvvj8Je9Xq+SXb77uWSU5C9bR8vuEuXEMOWqFlNwubi4uMyyx+MJKVfVfiKmyo/JJV8cyW4jq1TZKNnt+2kFygopV8WYqmM/EVPpmFwyFVrLq2JMx1s/RbNGlNVP4WKK9JibYBklWL5yosvIHaZMPxETMRFTdY2pPHElpAsXLtSsWbPUqVOnwLbOnTtrxowZeu+996KqY8+ePfJ4PEpLSwvZnpaWptzc3LD3yc3NDXt8cXGx9uzZE9iWl5enWrVqKSkpSYMGDdLTTz+tvn37RmzLlClTVKdOncAtPT1dkvTVV19JkrKyspSVlSVJ2rBhg7Zs2SJJWrdunbKzsyVJq1ev1q5duyRJK1euDHwlzvLlywNtW7p0qQ4cOCBJWrx4sQ4ePChJyszM1JEjR1RcXKzMzEwVFxfryJEjyszMlCQdPHhQixcvliQdOHBA47v6Ord9HaOxv/WVT61nNLqzr3x6I6PrOvgGY88mRle29ZX7NvNqSEtvlYxp6dKlknzjYvny5ZKknJycwJMTu3bt0urVqyX5Lt1et26dJN+nPW/YsIGYTtCYmteSJGnSaR6l1fCVHz3Do9QkKdntKye7pdQkX1mS0mr4jq+qMVXHfiKm0jE1r6UKreVVMabjrZ+iWSPK6qdwMUV6zB3c3KvBzX3lIS296tvMV76yrVc9mxj6iZiIiZiqZUybNm1SNCzjT51jULt2ba1YsUJdu3YN2b5u3Tr17t078ApjWXbv3q1mzZpp5cqVysjICGyfPHmyXnnlFW3cuLHUfdq3b6/hw4dr3LhxgW2ffPKJzjnnHOXk5KhJkyaSfFn7tm3b9PPPP+uDDz7QAw88oP/85z+lLuf1KygoUEFBQeD3/Px8paena9++fapXr14gu3e73SHl4uJiWZYVKLtcLrlcrojloqIiud3uQDkhIUGWZQXKku9ZjuByYmKijDGBstfrVaeJmSrwWnLJKMElFXotuSyjBMtXdltG7qCyy5KKgspbHhpc5WLyeDyBstfrVUJCQsSyx+ORMSZQDhcHMZ0YMXW4Z5G8spTsNir0+F7nOFr2JaUFHsmSlOSWCjyWLBkluaVNkwdXyZiqYz8RU+mYWt89X4kuxb2Wb35wQJWL6Xjrp1bjFpS5RhR4fH0TqZ82P3BhSExtJy4K+5hb5LUCr44WG0uJLiOvkTxB5a1TBtNPxERMxFTtYsrPz1f9+vWVl5en1NRURRLXe0jPO+883XbbbXrjjTfUtGlTSdL333+vMWPG6Pzzz4+qjoYNG8rtdpd6NfTHH38s9SqoX5MmTcIen5CQoAYNGgS2uVwutW3bVpLUtWtXZWVlacqUKRET0uTkZCUnJ5fa7na7Q36WLPs7M9pyYmJiXGXLsgJll8ulAq/v+2C9slToe5JVXmOp8NenFjzGkqecclWLyeVyRV2O1HZiOvFi8so3Fwo8R78jObTs+2lCylagXBVjiqdc1fuJmErH5JWlAv/6HcdaXlUfn46vfrLKXSPK6if/3wqOKVKfFZuj61KRN3yZfiImYiKm6hZT8O9lieuS3WeeeUYHDx5Uy5Yt1aZNG7Vt21atWrXSwYMH9fTTT0dVR1JSkrp3764lS5aEbF+yZInOOuussPfJyMgodfzixYvVo0ePkH9uScaYkFdAAQAAAADOi+sV0vT0dH3++edasmSJNm7cKGOMOnfurAsuuCCmesaOHaurrrpKPXr0UEZGhv75z39q586dge8VHTdunL7//nu9/PLLkqRRo0bpmWee0dixY3X99ddr1apVev755/XGG28E6pwyZYp69OihNm3aqLCwUJmZmXr55Zc1a9aseEIFAAAAAFSSmBPS4uJipaSkaP369erbt2+ZHxZUnmHDhmnv3r2aNGmScnJy1KVLF2VmZqpFixaSfG+e3blzZ+D4Vq1aKTMzU2PGjNGMGTPUtGlTTZ8+PfCVL5J06NAhjR49Wt99951q1Kihjh076tVXX9WwYcPibicAAC3vXhD3fbc/PMjGlgAAUH3EnJAmJCSoRYsWUX+Mb3lGjx6t0aNHh9334osvltrWu3dvff755xHre/DBB/Xggw/a0jYgkoqcmEqcnKL6YC4AAICKiOs9pBMnTtS4ceO0b98+u9sDAAAAADhBxPUe0unTp+vbb79V06ZN1aJFC5100kkh+8t6BRMAAAAAACnOhPSSSy6RZVmK4ytMAQAAAACQFGNCevjwYd155536z3/+o6KiIp1//vl6+umn1bBhw8pqHwAAKAPv4wUAHM9ieg/pfffdpxdffFGDBg3SFVdcoffff1833nhjZbUNAAAAAFCNxfQK6dy5c/X888/r8ssvlyRdeeWVOvvss+XxeOR2uyulgQAAAACA6immhHTXrl3q2bNn4PczzjhDCQkJ2r17t9LT021vHAAAOLa4BBgAcCzFdMmux+NRUlJSyLaEhAQVFxfb2igAAAAAQPUX0yukxhhde+21Sk5ODmw7cuSIRo0aFfLVL3PnzrWvhQAAAACAaimmhPSaa64pte0vf/mLbY0BAAAAAJw4YkpIX3jhhcpqBwAAAADgBBNTQgoAAIATFx96BcBuJKQAgCqjIie7nOgCAHD8ielTdgEAAAAAsAsJKQAAAADAESSkAAAAAABH8B5SAABwXOA9xgBQ/fAKKQAAAADAEbxCCgDiqwwAVB2sRwBOJLxCCgAAAABwBK+QAgAAVGO89zY+/N+AY4NXSAEAAAAAjiAhBQAAAAA4gkt2AQAAAMAhJ/oHmfEKKQAAAADAEbxCCgAAAEecKB8cdKK/AgaUhYQUACrBiXKSBQAAUBEkpACAaolXJAAAqPp4DykAAAAAwBEkpAAAAAAAR5CQAgAAAAAcQUIKAAAAAHAECSkAAAAAwBF8ym41VKPwiCTpl8RkybIkSYmeIiV4PPK43KEHHzr0651qSK5fn58oKpIKCyW3W0pJie/Yw4clY3zb3L/+zeJiqaDAd98aNeI79pdfJK9XSk6WEn4dvh6PdORIbMdallSz5tFjjxzx7UtKkhITyz22RuERFbkTVOz21WsZr1KKCn1/Nuno/yGpuEhur0fFbreK3ImBYwP/y5NOOlpvQYEv7sREXzsk3//l8GFfuWbNQH+qsND3v4/l2IQE3//CL1x/xnJsLH1vxzjx92fQsW6vR0nFRfJalgoSj7Y3uahALmNUmJAYGPMur0fJxUUylnQkMaVUvW6vp9xjk4sL5fJ6w/f9oUMh/VnmsSp/nIT0Z0XHSaT+jHecBPdnEP/acyQxScbyHZvgKVaip1hel0sFCUmBY1OKjsgy8s2zoL6vUXikVH/6jy1ISJTXVXbfh8z7X/n7M9ZxEiLMvI9lnMjr9bVNCunPSGtEeeMkwJgy1/ySY6pGUUGpYyP1fbj+DH4sKUxIDHtsQCWtEeXN+0jjJMQvv6hG4ZGo1ohI40SHDpV6fIhl3gcfG/DrOKlReCSqNSLQnyWO1aFDpdaIcH1fXn/+Ejy3fh0nScVFEfu+vHlfo/CIr21B60l5a0RU876SziNieSwJO++DzyOCxPP4oMLCo/0ZYT3hPCL284i4xknJ88KKnm96vUfnXIQ1othdvVM2XiGthrKmDlXW1KGq/0t+YNsN/5urrKlDdf+SWaEHN24s1aol7dx5dNuMGb5tI0aEHtuypW97VtbRbS++6Nt2+eWhx3bu7Nv++edHt735pm/bRReFHnv66b7tK1Yc3TZ/vm/bBReEHturl2/7okVHty1d6tuWkRF67IABvu3vvHN026ef+rb97nehxw4Z4tv+2mtHt335pW9bu3ahx151lbKmDtUVXywMbGqxP0dZU4fqfzOvCTn0oUXPKGvqUA1fMy+wrfHP+3z11q0bWu/Ysb7tDz10dFtenm9brVq+RdNvwgTftgkTjm4rLj56bF5eUCMe8m0bOzb079Wt69uek3N021NP+bb99a+hxzZr5tu+devRbf/8p2/bVVeFHtuunW/7l18e3fbaa75tQ4aEHvu73/m2f/rp0W3vvOPbNmBA6LEZGb7tS5cGNvXK/lxZU4fq/71+d8ihr705UVlTh+r8b1cHtp3x3dfKmjpU814q8X+46CKpVi0Nzloe2NTlh63KmjpU7z83OuTQp+c9qqypQzX0y/cD29ru2aWsqUN98yPII+89paypQ3X150e/eqRZ/k/KmjpUa5+5MuTY+5fMUtbUobrhf3OPbtyz52h/Bvv7333b7r//6LbDh48e6z+ZkHzH1Krlu08w/7F79hzd9thjvm033xx6bAxrxMezr1PW1KFqu2dXYNvQL99X1tShenreoyHHvv/caN//rcQakTV1qJ57+4GQY+e9NFZZU4fqjO++Dmw7/9vVypo6VK+9OTG0vWHWiLN2bFDW1KF659U7Qg596a37lDV1qPpvXhXY1m33JmVNHar3XrgltN4hQ5Q1dagu+eajwKaOP+1Q1tSh+uifN4Qc+uT8J0qtEdq61deuZs1Cjo20RmRNHaoNTw0LOXbi0ueUNXWoblr1/45uzMsLrPkJXk9g8x3LX1HW1KG6Y/krgW0JXk/g2NSCQ0GNCL9GbHhqmLKmDvWtWb8avmaesqYO1UOLngk59n8zr1HW1KFqsT9oPamkNaL/5lXKmjpUL711X8ih77x6h7KmDtVZOzYEtkVaI3TBBVGvEc+9/YCypg4ttUaoVi3fY12QstaIj2dfF3JsuDVCO3dKtWpFtUbU/yU/0J/B7v7ohVJrRI2igsCx/sRUkm775A1lTR3qu0+QcOcR/jWi5HnE2meuVNbUoWqW/1Ng29WfL1DW1KF65L2nQo79ePZ1pc4jylsjuvxw9HFncNbysGtEZZ1HxLJGzHrnoVJrRKTziHBrRHnnEXoq6H+Zk2PreUTLuxeo5d0L1PaueYFjTx3z/wLbp513rVSrll7OGBLY1vLuBSpKrXPcn0do0SLftl69Qo+94ALf9vnzj25bscK37fTTQ4/99TxCb755dNvnn4ddI3T55b7tL754dFtWVvRrRDVFQgoAAAAAcIRljDFON6Kqyc/PV506dZSXl6fU1FSnm1NKeV/2Xt4lu5sfv+TowVyy6xPjJbudxi+o0CW72fee5zuAS20qdKlNm7vmxX3J7vaHB4XU2+b+Dyp0yW7WAxeG9GeHO96J+5Ld7Q8POm4u2W35jw8Ch8ZzyW7WY5eE9H2nu/4b9yW72x8eFDLvW05cFNKfsY6TjU8EPRMfZt7HMk62PzQg5BI7/zoe7yW7Wx67xLfRGHX6m+9Vs3gv2d0+qW9I3/vbFu8lu9mP/MG3sRLWiJZ3L4j7kt1NT/zxaL2//KJOEzIrdMnutvvOD3l8aHn3grgv2f32sYt9G3+9FLPTPQsrdMnu5vv7hqwRLf8+P+5Ldrc/Mti38dc1ov19Syp0yW7WAxeGrCdt7/xv3JfsBtbxSjqP6PC3uXFfsrv94UEh5xEtJ30Y8Vip/Hmf9fAfKu2S3Zb3LgkcG884yXr80uP2PKKqXLLb6Q7f1XzxXrIbmAtVTLQ5VfW+IPkEFTyY/YrciUcfxIIFL2R+iYlHJ0m8xwZPPr+EhKMTNd5jgxcLP7c7fNtiOTal9P+srGNL/o+N5Qr7f/ct3Imljg1bb3Jy6OIt+RaycMcmJR19sInnWKnix8bS93aMkzD96XG59UuSu9T2kPcW/cob4Vh/vcHvr450bPDJkl+g70v8j8o8toRw4yRif1Z0nEiV1vfhYiuO8EAaOJELfj9kQkLYOkLe8/urSH0fbpxE7M9YxkmYeR/LOJEr/LyPtEbEMk6iXvMjHBup72N5LAlbbyWtEbHM+7LGSdT9GWGchOvPWOZ9WeOk5PGx9H1hQmLptsUyThShP38dJ8EJSaRjI837cGtluWtEkIj9WUnnEbGsEWH7M1K98Tw+BM/RCOuJHecRcY0TV9AFl8fZeUSVON90xbBGVFMkpABwginvKouyVNVnYQEAwPGJhBQAAACALXjSE7HiQ40AAAAAAI4gIQUAAAAAOIKEFAAAAADgCBJSAAAAAIAjSEgBAAAAAI4gIQUAAAAAOIKvfQGAKq4iH6Ev8TH6AFDd8NUqqE54hRQAAAAA4AheIcUJg2cTAQAAgKqFhBRAmUjkAQAAUFm4ZBcAAAAA4AheIQVwzPDhPAAAAAjGK6QAAAAAAEfwCilQBfA+TQAAAJyIeIUUAAAAAOAIElIAAAAAgCNISAEAAAAAjiAhBQAAAAA4goQUAAAAAOAIElIAAAAAgCNISAEAAAAAjuB7SAEAAAAgBhX5DnmJ75EPRkKKKouJDgAAALtU5NyS88rKwyW7AAAAAABHkJACAAAAABxBQgoAAAAAcATvIQXX0wMAAABwBK+QAgAAAAAcQUIKAAAAAHAECSkAAAAAwBEkpAAAAAAAR5CQAgAAAAAc4XhCOnPmTLVq1UopKSnq3r27VqxYUebxy5YtU/fu3ZWSkqLWrVtr9uzZIfufffZZ9ezZU/Xq1VO9evV0wQUXaPXq1ZUZAgAAAAAgDo4mpG+++aZuv/12TZgwQevWrVPPnj01YMAA7dy5M+zx2dnZGjhwoHr27Kl169Zp/PjxuvXWW/X2228Hjvnoo490xRVX6MMPP9SqVavUvHlz9evXT99///2xCgsAAAAAEAVHE9Inn3xSI0aM0MiRI9WpUydNmzZN6enpmjVrVtjjZ8+erebNm2vatGnq1KmTRo4cqeuuu06PP/544JjXXntNo0ePVteuXdWxY0c9++yz8nq9+uCDD45VWAAAAACAKDiWkBYWFmrt2rXq169fyPZ+/fpp5cqVYe+zatWqUsf3799fa9asUVFRUdj7HD58WEVFRapfv37EthQUFCg/Pz/kJkkejyfwM1y5uLg4pOz1esssFxUVhZSNMSFlY0ypsqSQstfrVbLLdz+XjJL8Zeto2V2inBimHBxHosvIbZlS5aQSZVdwWUfbbkdMweXi4uK4YkqwjBKC4gjXT9HElBxcdhtZKh1TsttI8u3zlRVSdsmU20/BsUYbU7h+quyxF0tMJfvJH1/weIsmpkj9VDKmcP10tFx+PwWPvYrMpwTLhF0jYp1PwWMvuJ9iiSlSP5Ucb7HOp+DxFjy34p1P/nLJsRfvfCo53sLNrVjHXsn5FG1MkfopeG7FM5+C+yl4vMU7nyKt5fHOp+CYgsdbpHJwP0V6nC35mBvvfEp2m1LrXkXmkytovFXGWh7vfPLHVHItj2c++csl+yne+RRpLY93Plkq/zyiomMv3vmU7DJhH3MrspZHGm92rOUVOTcqOd6Mse/cKJ41orzz8njmk7+fwp3DVmQtt+vcKKnEeKus8/KKzKfyOJaQ7tmzRx6PR2lpaSHb09LSlJubG/Y+ubm5YY8vLi7Wnj17wt7n7rvvVrNmzXTBBRdEbMuUKVNUp06dwC09PV2S9NVXX0mSsrKylJWVJUnasGGDtmzZIklat26dsrOzJUmrV6/Wrl27JEkrV65UTk6OJGn58uWBti1dulQHDhyQJC1evFgHDx6UJGVmZurIkSMqLi5WZmamiouLdeTIEWVmZkqSDh48qMWLF0uSDhw4oPFdfZ3bvo7R2N/6yqfWMxrd2Vc+vZHRdR18g7FnE6Mr2/rKfZt5NaSlt1RMQ1p61beZb/uVbb3q2cQ3eK/r4NXpjXzl0Z09OrWerzz2tx61r2NsjWnp0qWSfONi+fLlccU0uLlXg5t7AzGF66doYhrf1aPmtXxjY9JpHqXVUKmYHj3Do2S3lJokPXqGr41pNXzHS1LzWiq3n3bt2hV4f3O0MYXrp8oee7HEVLKf1q1bJ0nasmWLNmzYEHVMkfqpZEzh+unRMzxKTZKS3Sq3n4LHXkXm0+Dm3rBrRKzzKXjsBfdTLDFF6qecnJzAk33xzCd/P5Vc9+KdT/6YSo69eOfT6Y1MYD5lZ2cHxl6s8ym4n0rOp2hjitRPwetePPMpuJ/88ykrKyvu+RRpLY93PgXHFLyWB4+94HUvuJ+C14iyHnPjnU+TTvOUWvcqMp9OrWfCxmTXWh7vfPLHVHItj2c++cdeyX6Kdz5FWsvjnU9pNVTueURFx16882l8V0/Yx9yKrOXhYrJrLa/IudGjZ3hKnUfYdW4UzxpR3nl5PPPJ30/hzmErspbbdW409rce23MNu+bTpk2bFA3L+FPnY2z37t1q1qyZVq5cqYyMjMD2yZMn65VXXtHGjRtL3ad9+/YaPny4xo0bF9j2ySef6JxzzlFOTo6aNGkScvyjjz6qhx9+WB999JFOPfXUiG0pKChQQUFB4Pf8/Hylp6dr3759qlevXiC7d7vdIeXi4mJZlhUou1wuuVyuiOWioiK53e5AOSEhQZZlBcqS7+QsuJyYmChjTKDs9XrVaWKmCryWXDJKcEmFXksuyyjB8pXdlpE7qOyypKKg8paHBofE0W78fHmN5DGWEl0mUE5yGXmCysVG8vrLXskrS1se6GdLTB6PJ1D2er1KSEhQ63Hzo46pyGsFnlkq/jWOjQ8MCOmnNhMWRhVTssuoyF92GxV6pOyHB4fE1HHiAhV4JEtSklsq8FiyZAJll4wSXYrYT5sfGhwSa5tx86OKKVw/fftg/0ode75Yy48pXD9teuBCJSQkBMZbmwkLo4opUj8Fx1pUVKQO9ywK6ScTUvY9mETqp02TB4eMvdZ3z49rPvn7adODA0LWiDYTFsY8n/xjb1vQeLMsSx0mzI8qpkj9FG68xTqfgsdb8LrX4Z6FMc+n4H7aNmVQyNjzj7dY5lO48WaMUduJi2KeT8H9VHJutZ+4MKqYwvXTpsmDQ9a91nfPj3k+BffT5qDx1mHie3HNJ38/+cebf404Ot6in0/B/bSpxHiLVPb3U/AaUfJxtuRjbvB4i3Ut3/jgoJB1r9W4BXHNJ38/bf51vPljKjneYhl7W6cMDlmz24+fH9d88vfTxgcHhazlre6eH/N88o+9rZMvDOmnjve8F9d88vfTthKxdhg/P+b55O+nzQ/0L/M8oqJjr/3E9+I6N0p2GWU9ODDkMbfVuMy45pO/nzY/GH68xbOW+8ebfy1vP35+XOdGJmhuBZ9HtB63IK5zo0KvpW0PDajQGlHeeXn78fNjnk/+fvKPN/9a3u6exXGdG/n7ZuvkCwNtbztuftznRkkuo41B462yzsvjmU/5+fmqX7++8vLylJqaqkgSIu6pZA0bNpTb7S71auiPP/5Y6lVQvyZNmoQ9PiEhQQ0aNAjZ/vjjj+uhhx7S+++/X2YyKknJyclKTk4utd3tdof8LFn2d2a05cTExLjKlmUFyi6XSwVeS5JvAhX6nkCR11gq/PWpBY+x5CmnHBxH0a/1lSwXRlG2KyaXy1WqHGtMxSY0Dn+MwX0QTUwFwWVP+Fj92418C5WvbAXKXlkq8B4th+un4FijjSlcubLHnj/WaGIqGYe/DcHjLZqYIvVNyZi88u0L7qfQsu9npH4KGW8VmE/Fxgq7XsQ6nwoizK1YYorUT9GMt2jGXskx5vn1uFjnk78cvBYEb491PkUab7HOp0jjLSEh4eh4KycmX1m/xnG0HBxrPPMpuBw83vxx2bWWxzufgmOKtK4HlyM9tpb1mFveeCtr7FmWVSJWK+755A0ab5Wxlsc7n/RrHCVjjWc+hTtfcLvdgTbbtZbHO5+Myj+PqOjY84SZW9GOPf/fsmstj2W8xbqWV+TcKNzcMiFzK/5zo3jWiPLK8cwnXxzhz2ErspbbdW5UGGG82X1eHs98Cv69LI5dspuUlKTu3btryZIlIduXLFmis846K+x9MjIySh2/ePFi9ejRI+Sf+9hjj+mBBx7QwoUL1aNHD/sbDwAAAACoMEc/ZXfs2LF67rnnNGfOHGVlZWnMmDHauXOnRo0aJUkaN26crr766sDxo0aN0o4dOzR27FhlZWVpzpw5ev7553XHHXcEjnn00Uc1ceJEzZkzRy1btlRubq5yc3P1888/H/P4AAAAAACROXbJriQNGzZMe/fu1aRJk5STk6MuXbooMzNTLVq0kOR782zwd5K2atVKmZmZGjNmjGbMmKGmTZtq+vTpGjJkSOCYmTNnqrCwUEOHDg35W/fdd5/+8Y9/HJO4AAAAAADlczQhlaTRo0dr9OjRYfe9+OKLpbb17t1bn3/+ecT6tm/fblPLAAAAAACVyfGEFNVLy7sXVOj+2x8eZFNLAAAAAFR1jr6HFAAAAABw4iIhBQAAAAA4goQUAAAAAOAIElIAAAAAgCNISAEAAAAAjiAhBQAAAAA4goQUAAAAAOAIElIAAAAAgCNISAEAAAAAjiAhBQAAAAA4goQUAAAAAOAIElIAAAAAgCNISAEAAAAAjiAhBQAAAAA4goQUAAAAAOAIElIAAAAAgCNISAEAAAAAjiAhBQAAAAA4goQUAAAAAOAIElIAAAAAgCNISAEAAAAAjiAhBQAAAAA4goQUAAAAAOAIElIAAAAAgCNISAEAAAAAjiAhBQAAAAA4goQUAAAAAOAIElIAAAAAgCNISAEAAAAAjiAhBQAAAAA4goQUAAAAAOAIElIAAAAAgCNISAEAAAAAjiAhBQAAAAA4goQUAAAAAOAIElIAAAAAgCNISAEAAAAAjiAhBQAAAAA4goQUAAAAAOAIElIAAAAAgCNISAEAAAAAjiAhBQAAAAA4goQUAAAAAOAIElIAAAAAgCNISAEAAAAAjiAhBQAAAAA4goQUAAAAAOAIElIAAAAAgCNISAEAAAAAjiAhBQAAAAA4goQUAAAAAOAIElIAAAAAgCNISAEAAAAAjiAhBQAAAAA4goQUAAAAAOAIElIAAAAAgCNISAEAAAAAjiAhBQAAAAA4goQUAAAAAOAIElIAAAAAgCNISAEAAAAAjiAhBQAAAAA4goQUAAAAAOAIElIAAAAAgCNISAEAAAAAjiAhBQAAAAA4goQUAAAAAOAIElIAAAAAgCNISAEAAAAAjiAhBQAAAAA4wvGEdObMmWrVqpVSUlLUvXt3rVixoszjly1bpu7duyslJUWtW7fW7NmzQ/Z//fXXGjJkiFq2bCnLsjRt2rRKbD0AAAAAIF4JTv7xN998U7fffrtmzpyps88+W//3f/+nAQMG6JtvvlHz5s1LHZ+dna2BAwfq+uuv16uvvqpPPvlEo0ePVqNGjTRkyBBJ0uHDh9W6dWv96U9/0pgxY451SAAA4Fct715Qoftvf3iQTS0BAFRVjr5C+uSTT2rEiBEaOXKkOnXqpGnTpik9PV2zZs0Ke/zs2bPVvHlzTZs2TZ06ddLIkSN13XXX6fHHHw8cc/rpp+uxxx7T5ZdfruTk5GMVCgAAAAAgRo4lpIWFhVq7dq369esXsr1fv35auXJl2PusWrWq1PH9+/fXmjVrVFRUFHdbCgoKlJ+fH3KTJI/HE/gZrlxcXBxS9nq9ZZaLiopCysaYkLIxplRZUkjZ6/Uq2eW7n0tGSf6ydbTsLlFODFMOjiPRZeS2TKlyUomyK7iso20Pjsn6dXuy2/xaNkp2+35agbJCyi6ZQExer1fFxcVxxZRgGSUExRGun6KJKTm4HIgjtJ9iiSlSPwXHGm1M4fqpssdeLDGV7Cd/fMHjLZqYIvVTyZjC9VMsY8/r9QZirch8SrBM2DUi1vkUPPaC+yne+RQcU8nxFut8Ch5vwXMr3vnkL5cce/HOp5LjLdzcinXslZxP0cYUqZ+C51Y88ym4n4LHW7zzKdJaHu98qoy1vORjbrzzKdltSq17FZlPrqDxVhlrebzzyR9TybU8nvnkLwfPJ4/HE/d8irSWxzufLIWeGwWXw/VNcLlkTJHO9+KdT8kuE/YxtyJreaTxZsdaXpFzo5LjzRj7zo2i7adYzsvjmU/+fgp3Xl6Rtdyuc6OkEuPNjlzD7vlUHscS0j179sjj8SgtLS1ke1pamnJzc8PeJzc3N+zxxcXF2rNnT9xtmTJliurUqRO4paenS5K++uorSVJWVpaysrIkSRs2bNCWLVskSevWrVN2drYkafXq1dq1a5ckaeXKlcrJyZEkLV++PNC2pUuX6sCBA5KkxYsX6+DBg5KkzMxMHTlyRMXFxcrMzFRxcbGOHDmizMxMSdLBgwe1ePFiSdKBAwc0vquvc9vXMRr7W1/51HpGozv7yqc3Mrqug28w9mxidGVbX7lvM6+GtPSWimlIS6/6NvNtv7KtVz2b+AbvdR28Or2Rrzy6s0en1vOVx/7Wo/Z1TNiY0mr4/qePnuFRapKU7PaVk91SapKvLElpNaRJp/nKzWspENOePXu0fPnyuGIa3Nyrwc29gZjC9VM0MY3v6lHzWr44Jp3mCcQU3E+xxBSpn3bt2qXVq1fHFFO4fqrssRdLTCX7ad26dZKkLVu2aMOGDVHHFKmfSsYUrp9iGXsHDhzQ0qVLY4opUj+FWyNinU/BYy+4n+KdT8Ex5eTkBJ7si2c++fup5LoX73zyx1Ry7MU7n05vZALzKTs7OzD2Yp1Pwf1Ucj5FG1Okfgpey+OZT8H95J9PWVlZcc+nSGt5vPOpMtbyko+58c6nSad5Sq17FZlPp9YzgflUGWt5vPPJH1PJtTye+eTvp+D5tGXLlrjnU6S1PN75lFZDIedG/rU8eOwFr3vB/VQypuD5FLyWxzufxnf1hH3MrchaHi4mu9byipwbPXqGp9R5hF3nRtH2Uyzn5fHMJ38/hTsvr8habte50djfemzPNeyaT5s2bVI0LONPnY+x3bt3q1mzZlq5cqUyMjIC2ydPnqxXXnlFGzduLHWf9u3ba/jw4Ro3blxg2yeffKJzzjlHOTk5atKkScjxLVu21O23367bb7+9zLYUFBSooKAg8Ht+fr7S09O1b98+1atXL5Ddu93ukHJxcbEsywqUXS6XXC5XxHJRUZHcbnegnJCQIMuyAmXJd3IWXE5MTJQxJlD2er3qNDFTBV5LLhkluKRCryWXZZRg+cpuy8gdVHZZUlFQectDg0PiaDd+vrxG8hhLiS4TKCe5jDxB5WIjef1lr+SVpS0P9AuJqf09i2RkKdltVOiRjHwTusAjWZKS3FKBx5IlEyi7ZJTokgq8lrY9NEBer1cJCQlqPW5+1DEVea3AM0vFv8ax8YEBIf3UZsLCqGJKdhkV+cu/xpH98OCQfuo4cUHUMYXrp80PDZbX6w3E2mbc/KhiCtdP3z7Yv1LHni/W8mMK10+bHrhQCQkJgfHWZsLCqGKK1E/BsRYVFanDPYtC+imWsbdpsq8PPB6PEhMT1fru+XHNJ38/bXpwQMga0WbCwpjnk3/sbQsab5ZlqcOE+XHNJ38/hRtvsc6n4PEWvO51uGdhzPMpuJ+2TRkUMvb84y2W+RRuvBlj1HbiopjnU3A/lZxb7ScujCqmcP20afLgkLW89d3zY55Pwf20OWi8dZj4Xlzzyd9P/vHmXyOOjrf41vJNQeOt7cRFcc0nfz9tnXxhyGNu8HiLdS3f+OCgkHWv1bgFcc0nfz9t/nW8Bcca63zyl7dOGRyyZrcfPz+u+eTvp40PDgpZy1vdPT/m+eTvp62TL5QxJjC3Ot7zXlzzyd9P20rE2mH8/Jjnk7+fNj/QP3Bu5F/Lg9e6SGX/GhH8+FTyHM/j8aj9xPfiOjdKdhllPTgw5DG31bjMuM+NXDLa/GD48RbPWu4fb/61vP34+XGdG5mguRV8HtF63IK4zo0KS5wHRttPsZyXtx8/P+b55O8n/3jzr+Xt7lkc17mRv2+C17e24+bHfW6U5DLaGDTe7Mg17JpP+fn5ql+/vvLy8pSamqpIHPtQo4YNG8rtdpd6NfTHH38s9SqoX5MmTcIen5CQoAYNGsTdluTk5LDvN3W73SE/S5b9nRltOTExMa6yZVmBssvlUoHXkuSbQIW+J1DkNZYKf31qwWMsecopB8dR9Gt9JcuFUZRLttfIt6/Ac/SYgl9frTchZStQ9spSwa9x+BeMeGIqNqFx+GMM7oNoYioILnvCx+rfHk1MkfopONZoYwpXruyx5481mphKxuFvQ/B4iyamSH1TMiZv2PEW/dgLGW8VmE/Fxgq7XsQ6nwoizK1451NwTNGMt2jGXskx5vn1uFjnk78cvL4Fb491PkUab7HOp0jjLSEh4eh4KycmX1m/xnG0HBxrPPMpuBw83vxx2bWWxzufKmMtL/mYW954K2vsWZZVIlYr7vnkDRpvlbGWxzuf9GscJWONZz6FO19wu92BNtu1lsc7n4xCz438fRBNOdJ5XcmyJ8zcinbs+f+WXWt5LOMt1rW8IudG4eaWCZlb8Z8bRdtP4WKKVI5nPvniCH9eXpG13K5zo8II460iuYZd8yn497I4dsluUlKSunfvriVLloRsX7Jkic4666yw98nIyCh1/OLFi9WjR4+Qfy4AAAAAoOpz9FN2x44dq+eee05z5sxRVlaWxowZo507d2rUqFGSpHHjxunqq68OHD9q1Cjt2LFDY8eOVVZWlubMmaPnn39ed9xxR+CYwsJCrV+/XuvXr1dhYaG+//57rV+/Xt9+++0xjw8AAAAAEJmj30M6bNgw7d27V5MmTVJOTo66dOmizMxMtWjRQpLvzbM7d+4MHN+qVStlZmZqzJgxmjFjhpo2barp06cHvoNU8r03tVu3boHfH3/8cT3++OPq3bu3Pvroo2MWGwAAAACgbI4mpJI0evRojR49Ouy+F198sdS23r176/PPP49YX8uWLeXQ5zQBAAAAAGLg6CW7AAAAAIATFwkpAAAAAMARJKQAAAAAAEc4/h5SAAAAACip5d0LKnT/7Q8PsqklqEy8QgoAAAAAcASvkALVDM8mAgAA4HjBK6QAAAAAAEeQkAIAAAAAHEFCCgAAAABwBAkpAAAAAMARJKQAAAAAAEeQkAIAAAAAHEFCCgAAAABwBAkpAAAAAMARJKQAAAAAAEeQkAIAAAAAHEFCCgAAAABwBAkpAAAAAMARJKQAAAAAAEeQkAIAAAAAHEFCCgAAAABwBAkpAAAAAMARJKQAAAAAAEeQkAIAAAAAHEFCCgAAAABwBAkpAAAAAMARJKQAAAAAAEeQkAIAAAAAHEFCCgAAAABwBAkpAAAAAMARJKQAAAAAAEeQkAIAAAAAHEFCCgAAAABwBAkpAAAAAMARJKQAAAAAAEeQkAIAAAAAHEFCCgAAAABwBAkpAAAAAMARJKQAAAAAAEeQkAIAAAAAHEFCCgAAAABwBAkpAAAAAMARJKQAAAAAAEeQkAIAAAAAHEFCCgAAAABwBAkpAAAAAMARJKQAAAAAAEeQkAIAAAAAHEFCCgAAAABwBAkpAAAAAMARJKQAAAAAAEeQkAIAAAAAHEFCCgAAAABwBAkpAAAAAMARJKQAAAAAAEeQkAIAAAAAHEFCCgAAAABwBAkpAAAAAMARJKQAAAAAAEeQkAIAAAAAHEFCCgAAAABwBAkpAAAAAMARJKQAAAAAAEeQkAIAAAAAHEFCCgAAAABwBAkpAAAAAMARJKQAAAAAAEeQkAIAAAAAHEFCCgAAAABwBAkpAAAAAMARJKQAAAAAAEc4npDOnDlTrVq1UkpKirp3764VK1aUefyyZcvUvXt3paSkqHXr1po9e3apY95++2117txZycnJ6ty5s955553Kaj4AAAAAIE6OJqRvvvmmbr/9dk2YMEHr1q1Tz549NWDAAO3cuTPs8dnZ2Ro4cKB69uypdevWafz48br11lv19ttvB45ZtWqVhg0bpquuukpffPGFrrrqKl122WX63//+d6zCAgAAAABEwdGE9Mknn9SIESM0cuRIderUSdOmTVN6erpmzZoV9vjZs2erefPmmjZtmjp16qSRI0fquuuu0+OPPx44Ztq0aerbt6/GjRunjh07aty4cTr//PM1bdq0YxQVAAAAACAajiWkhYWFWrt2rfr16xeyvV+/flq5cmXY+6xatarU8f3799eaNWtUVFRU5jGR6gQAAAAAOCPBqT+8Z88eeTwepaWlhWxPS0tTbm5u2Pvk5uaGPb64uFh79uzRySefHPGYSHVKUkFBgQoKCgK/5+XlSZL2798vSfJ4PJIkt9sdUi4uLpZlWYGyy+WSy+WKWC4qKpLb7Q6UExISZFlWoCxJxcXFIeXExEQZYwJlr9erhKJDKvRacsnI7ZKKvJZclpHb8pXdlpGrjHJ+fn5IHO6iQ/IayWMsJbpMxLLHSF5/2St5ZWnv3r0hMZmCQzKylOQ2KvJIRlKSWyr0SJakRLdU6LFkyQTKLhkluKRCr6UDBw74YkxIkAoPRR1TkddSgmV8/7df27h///6QfvIWHI4qpiSXUbG//Gsc+fn5If2UUHwo6pjC9VN+fr68Xm8gVqvwUFQxheubffv2hYw3FR6KKqZI/ZSXlxcy9nyxlh9TuL7Zt2+fEhISAuPNW3A4qpgi9VNwrEVFRVLBoahiCtdP/j7weDxKTEyUCg7FNZ/8/eQfb8Gxxjqf/P0UPN4sy1JC8aG45pO/n6Idb9GMvX379oWse1bhoZjnU3A/+cdb6NyKbT6FG2/GGHkLDsc8n8qcW/7xFsN88vdTfn5+yFqugkMxz6fgcvB4cxUeims++fvJP978a/nR8RbfWh483rwFh+OaT/5+8j8Ohxtvsa7leXl5IY+53oJDcc0nfz/5x1twrLHOJ385Pz8/5HzBXXQorvnk7yf/eUxgLS84FPN88rd9//79vseCX+eWu+hQXPPJ308lY00oOhTzfPL30969ewPnRv61PHiti1T2rxHBj08lz/E8Hk/E8RbN2Dtw4EDI+Z634HDc50YuxTbeyht7/j7wr+UlzwNjHXvB4y0hISFkvMW6lgefBwaPt3jX8rBzK45zI0sKjDf/Wu4tOBzXuVHw3PKPN1fhobjPjRJLjDc7cg275lN+fr4kyRijMhmHfP/990aSWblyZcj2Bx980HTo0CHsfdq1a2ceeuihkG0ff/yxkWRycnKMMcYkJiaa119/PeSYV1991SQnJ0dsy3333WfkG3fcuHHjxo0bN27cuHHjxs2m265du8rMCx17hbRhw4Zyu92lXrn88ccfS73C6dekSZOwxyckJKhBgwZlHhOpTkkaN26cxo4dG/jd6/Vq3759atCggSzLiikup+Xn5ys9PV27du1SampqlamLtlWN+mhb9aqLtlWN+mib83XRtqpR34nSthMlTtpWNeqzu23HkjFGBw8eVNOmTcs8zrGENCkpSd27d9eSJUt06aWXBrYvWbJEF198cdj7ZGRk6N133w3ZtnjxYvXo0cN3CdSvxyxZskRjxowJOeass86K2Jbk5GQlJyeHbKtbt26sIVUpqamptg1aO+uyuz7a5nxddtd3orTtRInT7vpom/N12V0fbXO+LrvrO1HadqLEaXd9tM35uo6lOnXqlHuMYwmpJI0dO1ZXXXWVevTooYyMDP3zn//Uzp07NWrUKEm+Vy6///57vfzyy5KkUaNG6ZlnntHYsWN1/fXXa9WqVXr++ef1xhtvBOq87bbb1KtXLz3yyCO6+OKL9d///lfvv/++Pv74Y0diBAAAAACE52hCOmzYMO3du1eTJk1STk6OunTposzMTLVo0UKSlJOTE/KdpK1atVJmZqbGjBmjGTNmqGnTppo+fbqGDBkSOOass87Sv/71L02cOFH33HOP2rRpozfffFNnnnnmMY8PAAAAABCZowmpJI0ePVqjR48Ou+/FF18sta137976/PPPy6xz6NChGjp0qB3NO+4kJyfrvvvuK3UJstN12V0fbXO+LrvrO1HadqLEaXd9tM35uuyuj7Y5X5fd9Z0obTtR4rS7PtrmfF1VlWVMeZ/DCwAAAACA/VxONwAAAAAAcGIiIQUAAAAAOIKEFAAAAADgCBJSAAAAAIAjSEiBIIcOHdLy5cudbgZskpWVpdatWzvdDAAAAETg+Ne+oOLeeustvfHGG9q8ebMsy1K7du305z//2bavvikuLtbu3bvVvHlzW+qrrDrt8O2336pPnz7yeDxONwU2KCws1I4dO2ytc//+/Xr33Xd19dVX21pvLLxer1yu0s8ner1efffdd1HPK2OMtm/frvT0dCUkJKiwsFDvvPOOCgoKNHDgQDVs2NDupgNApZs3b165xyQkJKhJkybq0qWLkpKS4v5b8T4m2LWOh/PDDz+ooKAgrjp4XIhNUVGRJkyYoLlz56p+/fq68cYbNXz48MD+H374QU2bNuW8sjwGxy2Px2Muu+wyY1mW6dChg7n44ovNRRddZNq3b29cLpcZNmyY8Xq9Ff4769evNy6Xy4YWx1fnjBkzzPnnn2/+9Kc/mQ8++CBk308//WRatWrlSLuitW/fPvPSSy/ZUtfPP/9sli1bZktd3377renTp48tdcWruLg45PdPP/3ULFu2zBQWFtpSv9Nj1255eXnmT3/6k0lJSTGNGzc29957b8j/MDc3N+q2bdy40bRo0cK4XC7Ttm1bs23bNtO9e3dz0kknmZo1a5qGDRuazZs3x9xGj8cTcfuOHTtiri+SWOdVZawjdsXq9XrNtm3bTFFRkTHGmIKCAvOvf/3LvPTSS+ann36KqU2HDh0yo0ePNk2bNjWNGjUyV1xxRcx1RCvWPti9e7d55ZVXzIIFC0xBQUHIvp9//tncf//9MbfBzvFmZz/YWRdiZ1lW1LeTTz7ZLF++PO6/Fetjgp3reH5+vrnyyitN8+bNzdVXX20KCgrM6NGjjWVZxuVymV69epm8vLyo21ZZjwsl9enTx2zfvj3m+y1evNjce++9gTV82bJl5sILLzR9+vQxc+bMqXC74nHfffeZtLQ089hjj5kJEyaYOnXqmBtuuCGwPzc311iW5UjbjickpMexJ554wtSvX9+8++67pfb997//NfXr1zdTp06t8N9x8qT+qaeeMjVr1jQ33XST+ctf/mKSk5PNQw89FNgfy8JtjDH16tUr85aamlqlExgn6yosLDR33nmnadOmjTn99NNLLf6x9MXu3bvN2Wefbdxut+nVq5fZt2+fGTRoUOAEoX379mb37t0xxRNOPP+vvLy8Mm8rVqywdYx88803USdDt956q2nfvr156623zLPPPmtatGhhBg0aFDi5j+WBz/8E1oYNG8ztt99uOnfubC6++GJTWFhoCgoKzMUXX2z+8pe/RB2HnSdZ0Yilb+1eR6ryEwN33HGHqVmzprn++uvNrbfeaho2bGiGDh0a9f1jEUsfrF692tStW9ekpqaaGjVqmHbt2pmvvvoqsN/JPjDG3n44Vif1xhizc+dOM3z48Jjuc/jwYbNixQrz9ddfl9r3yy+/xPQkw/r1680DDzxgZsyYUSrRzsvLi7lt33zzjZkzZ47JysoyxhiTlZVlRo0aZYYPH17qiaSK8nq9Jicnx9x0002mW7duEY+z+zHBznX85ptvNh07djTTp0835557rrn44otNly5dzMcff2yWL19uunTpYsaPHx912+x+XPjvf/8b9uZ2u80zzzwT+D0ar7zyiklISDCnnXaaqVWrlnnhhRdM3bp1zciRI82IESNMUlKSeeutt6JuW3miXd/atm0bch7+7bffmnbt2plrr73WeL3euB/75s+fb0aMGGHuvPPOwHzw27dvn+MvKNiNhPQ49tvf/tY8//zzEfc/99xzpkuXLuXW061btzJvHTt2jHky2VVn586dzWuvvRb4feXKlaZx48bmnnvuMcbEftJRs2ZN87e//c28+OKLYW/3339/lU5gYj0JL+t21113xdQuO58FvOqqq8xZZ51l5s2bZ4YNG2bOOuss07NnT/Pdd9+ZnTt3mp49e5qbbrop6rZFEk9C6n9mOdLNv98usbSxefPm5sMPPwz8vmfPHnPmmWeafv36mSNHjsQ0Hxo1amTWrVtnjPG9OmVZllmxYkVg/8qVK03z5s2jjsPOkyxj7J1Xdq8jVfmJgdatW5s33ngj8Pv//vc/k5CQUOpqhGjY2QcXXHCBue6664zH4zH5+flm9OjRpkGDBubzzz83xjjbB8bY2w9292lZYl3jNm3aZFq0aBFYx3r37h3y5F8s/bBo0SKTlJRkTjnlFNO8eXPTsGFDs3Tp0rjqMsaY9957zyQlJZn69eublJQU895775lGjRqZCy64wJx//vkmISHB9qTUGGOys7NNcnJyxP12PybYuY6np6cH/ufff/+9sSzLzJs3L7B/wYIFpkOHDlG3ze7HBf//pqxXqKONtWvXruapp54yxhjz/vvvmxo1apgnn3wysP+JJ54wZ599dtRtK8/69eujWkNq1KhhsrOzQ7Z9//33pkOHDubKK68033//fcznDK+99ppxu91m0KBB5pxzzjEpKSnm1VdfDey3+wneqsAyxhinLxtGfGrUqKFNmzZFfI/Ajh071LFjR/3yyy9l1pOSkqLLL79crVq1Crs/JydHzz77bEzXv9tVZ82aNfXNN9+oZcuWgW1ff/21zj//fA0fPly33357TNfmn3322brssst02223hd3/xRdf6LTTTospVpfLJcuyIu43xsiyrKjqrF+/fpn7PR6Pfv7556jqcrlcOvnkkyO+N6awsFC5ublRx9quXTtNnTpVgwcPliRt3bpVAwYM0Nlnn605c+boxx9/jLovmjZtqrlz5+r3v/+99u3bp4YNG2rJkiU6//zzJUkffvihRo4cqa1bt5ZZT7169cr83xcXF+vQoUMx9WedOnU0YcIEnXnmmWH3b9myRX/961+jrnPs2LFl7v/pp5/0+uuvR1XfSSedpK+++ipkXh08eFD9+/dXjRo19Nxzz6lt27ZR1VWzZk1t3LgxsH7Url1b69evV5s2bSRJu3btUrt27XTkyJFy65KkFi1a6KWXXtK5554rSdq7d68GDRqkOnXqaN68eTpw4EBMc9XOeWX3OmJnrI0bN9bixYvVtWtXHTp0SLVr19by5ct1zjnnSJJWrVqlyy+/POr3QiclJSk7O1vNmjULbKtRo4Y2b96s9PT0qOrws3tt+/TTT9W+ffvAtkcffVQPP/ywFi1apObNmzvWB5K9/WBnXeW9F3Lbtm3629/+FnWcl156qYqLi/XCCy/owIEDGjt2rL766it99NFHat68eUzvdzvrrLPUp08fTZ48WcYYPf7445o0aZLeeustXXjhhTG/d+6ss87SeeedpwcffFD/+te/NHr0aN14442aPHmyJGnChAn67LPPtHjx4qjqk6J/n2ZeXp7q1KkTtg67HxPsXMdTUlK0ZcuWwNw+6aSTtG7dusA827Fjhzp37qxDhw5F1Ta7HxcGDBggt9utOXPmqHHjxoHtiYmJ+uKLL9S5c+eo6pGkWrVq6csvvwz835KSkrRmzRqdeuqpkqRNmzbp7LPP1p49e6Kq749//GOZ+/Py8vTRRx+V2w+tW7fWs88+Gzh/8du9e7f69Omj5s2ba+nSpTGdh5x22mkaPny4brnlFknSv//9bw0fPlzTpk3TiBEjquX7UvlQo+NYjRo1dODAgYgJaX5+vmrUqFFuPV26dNGZZ56pG2+8Mez+9evX69lnn42pbXbV2bBhQ+3atSvkRPKUU07R0qVLdd555+n777+PqV2DBg3SgQMHIu6vX79+zB9MULt27agerKJRUFCgG2+8Ub/97W/D7t+xY4fuv//+qOpq0aKFHnnkEV122WVh969fv17du3ePqi5J+v7779WlS5fA723atNFHH32k8847T1dddZUeffTRqOvav39/4IS5fv36qlmzplq0aBFSd05OTrn1TJs2Leq/Ga3TTjtNktS7d++w++vWratYnsd76qmn1LVrV6Wmpobd//PPP0ddV3p6urKyskJOZGrXrq3FixerX79+uvTSS6Ouq2nTptq5c2dg/Xj00UdDThh++ukn1atXL+r69uzZE9KHDRo00JIlS9S/f38NHDhQzz33XNR1SfbOK7vXETtj/fnnnwNPRJ100kk66aSTdPLJJwf2/+Y3v9EPP/wQdX0ej6fUk1AJCQkqLi6Oug4/O/tAUqmT2Lvuuksul0v9+vXTnDlzYmqb3ePNzn6ws65LLrlElmWVueaU9aRBSStXrtT777+vhg0bqmHDhpo3b55uuukm9ezZUx9++KFOOumkqOv6+uuv9corrwTacOedd+o3v/mNhg4dqjfeeENnnHFG1HX563v55ZclSZdddpmuuuoqDRkyJLD/iiuu0PPPPx9VXfn5+Ro5cqTeffddpaamatSoUbr33nvldrsl+da3Vq1aBU7oIyWjkv2PCXau4w0aNNBPP/0USEgvvvhi1a1bN7D/559/VnJyctT12f248N5772nq1Kk6/fTTNWPGjMAT2vFITExUYWFh4Pfk5GTVqlUr8HtSUlK5L8AEe/fdd9W3b1+lpaWF3R9tsnfeeefp9ddfL5WQNm3aVEuXLg08aRaLzZs3h/yvhg4dqoYNG+qiiy5SUVFRTGPkeEFCehzLyMjQrFmzNGvWrLD7Z8yYoYyMjHLrOeecc7Rp06aI+2vXrq1evXrF1Da76jznnHP09ttvq2fPniHbO3furA8++EB9+vSJqV3jx48vc396erpeeOGFmOq088Gqa9euSk9P1zXXXBN2/xdffBF1Qtq9e3etXbs2YkJa3klOSU2aNNHWrVtDTur9C26fPn0itjmcxo0bKycnJ/AgevPNN4e8Orx///6oToxi+ZvR+vOf/6zDhw9H3N+kSRPdd999UdfXrl07jRkzRn/5y1/C7o/liYG+ffvqhRde0MCBA0O216pVS4sWLVLfvn2jbtcFF1ygjRs3Bl61Kfnk0eLFiwNjOxp2nmRJ9s4ru9eRqvzEgDFG1157bchJ6JEjRzRq1KiQOTV37txy67KzD7p06aKVK1cGXs3wu+OOO2SM0RVXXBFVPX52jzc7+8HOuk4++WTNmDFDl1xySdj9sT6x+MsvvyghIfTUb8aMGXK5XOrdu7def/31qOtKTk4u9QTvFVdcIZfLpcsvv1xPPPFE1HWV5HK5lJKSEpJc1a5dW3l5eVHd/5577tEXX3yhV155RQcOHNCDDz6otWvXau7cuYEnbKIdu3Y/JvTr18+2dfzUU0/VZ599FpirJfvvs88+U6dOnaKuz+7HBUkaM2aMzjvvPP35z3/Wu+++q6lTp8Z0f7+2bdtq48aN6tChgyTfk+S1a9cO7N+6dat+85vfRF1fp06dNGTIEI0YMSLs/vXr12v+/Pnl1nPPPfdo48aNYfc1a9ZMy5cvj+lVfUlKTU3VDz/8ELK+nXvuuXr33Xc1ePBgfffddzHVd1w49lcJwy6ffPKJSUxMNH/605/M//73P5OXl2cOHDhgVq1aZYYOHWoSExPNxx9/bMvf8r+n4Fj74osvyvzktK+//tpcfPHFUdc3fPhwk5+fb0PLjvrnP/8ZeF9DOLm5ueYf//hHVHVNnjy5zGN37txprr322qjq+vrrr81nn30WcX9hYWFMn3I3YsQIc91114Xd991335m2bdtG/Z6Giy66yEybNi3i/meeecacd955UbetKvvzn/9sbr/99oj7o32fijG+DzII90EkfgcPHjQfffRRzG0MJzs7O6YPlrrlllsifnhOfn6+OfPMM2N6z4ud86q8deSrr76Kui5j7I31r3/9q3n22Wcj7p8yZYoZOHBg1G275pprzLXXXlvuLRr//Oc/y5ynsfTBs88+W+b7Jh955BHTsmXLqOoyxvdhLnaONzv7wc66/vCHPwTe6xxOLOuHMcacfvrp5uWXXw6776abbjJ169aN+v/Wt29f89hjj4Xd9/rrr5vExMSY+uDUU0817733XuD3L7/8MvApxcYYs2LFiqg/AM7O92nabd++fSEf6FVSLOv43r17zf79+yPuz8zMDPk/VNS2bdvi/sDBw4cPm7/+9a+mXbt2xu12l/lYFs7cuXPL/JaBKVOmmIkTJ0Zd37XXXmtGjx4dcf8333wT1ZpUGeeVF198sbn33nvD7vvwww/NSSedVO3eQ0pCepybO3euadiwYak32Ddo0MD8+9//rlDdBw4cMDNmzDDdunWrcgM/3ra5XC7zww8/VGLLwvvuu++O+d+02/bt283ChQsj7t+9e7d58cUXbflbq1evNl9++WW5x9WtW7fcT06uV69eTH/7gw8+MJ06dQr7UfkHDhwwnTt3jukrAnJycuL6ePt42xbt1wLZHWc0J1nBH8pwPLPzhLI8sZ4Abt26NeJXoVQnx7IPjPH1Q05Ojm117dy5M6pjly9fHpKklfTzzz/HFOdDDz1kBgwYEHH/jTfeGHWCO3fu3DKfbHv99dfNueeeG3XbZs2aZebPnx9x//jx482IESOiqqtmzZpm27ZtIdvy8/NNRkaGOe+888y2bduiPnf45ZdfQj5F9e677zZjxowJ3O644w7zyy+/RFWXMb7ErKz67rzzzpjqs5PdsYbz3//+19x+++2BczGnzo+OHDliDh06VOF6KuO88qOPPgr5JPiSPvzww6ifWDxekJBWA4cOHTJz5841jzzyiHnkkUfMO++8U6FJ9sEHH5grr7zS1KhRw3Ts2NFMmDAh8EmI0dq8ebO5/PLLI57sXnHFFWbr1q3HvG2WZR3ThDQnJ8fcfPPNJiUlJab77dmzJ1DeuXOnueeee8wdd9xRoe9Ks4udbYtUV7RJVaRPSy55i8Uf/vCHkE/uK+mpp54yl1xySUx12sXOth3LOOOdBx6Pxzz//PNm0KBB5pRTTjFdunQxf/jDH8xLL70U03csV8Z6ZFfb7K6r5MnRZZddZnJzc2Oqw2/AgAHmwIEDgd8ffPDBkFdj9uzZYzp16hRX3RVld9uOVayxzoWtW7fa8n3iwfVV1Scs7Gxbhw4dzIIFC0ptP3jwoMnIyDC/+93vok5IZ8+ebQYPHhz4vVatWubMM8805557rjn33HNNkyZNylxLK7M+u9c2u2MtS7yPC3Yp71PE/bfyVMZ5pV1tO56QkB7H7HyWbdeuXeaBBx4wrVq1Mo0bNzY333yzSUhIiPmSCr/rr7/e3HnnnRH333XXXWbUqFHHvG2WZZkff/wx5vuVZf/+/ebPf/6zadiwoTn55JPNU089ZTwej7nnnntMjRo1TI8ePczrr78eVV0bNmwIfH9dhw4dzLp160xaWpqpVauWSU1NNW6327zzzjtR1dW1a9dyv36nrO9eq8y22VlXeYIv+YpG8+bNzTfffBNxf1ZWlklPT4+6vkhfGVC3bl1z5plnmrffftuRttkdp53zwBjfdwQOHDjQWJZlunbtai6//HIzbNgwc+qppxrLsmK6VN/O9cjfNv/35la0bXbGaUzpk6NatWrF9eRfuLpq164dUlcslz3avR6VTLwr0jZj7I3Vzrlg5xMMdtdn59pmd9vsvKy+Z8+eZu7cuYHfS86pV155xfz+97+Pum121mf32mZ3rHbOBbvXELu+zqcyzivLa5v/Vp3woUbHsZdfflnz588PfBLXM888o1NOOSXwybobN27UySefrDFjxpRZz8CBA/Xxxx9r8ODBevrpp3XhhRfK7XZr9uzZcbdt+fLlgU/fC+eyyy7Tn//853LrqYy2tW/fvtxPJdy3b1/U9Y0fP17Lly/XNddco4ULF2rMmDFauHChjhw5ovfeey/iB4KEc9ddd+m3v/2tXn31Vb366qsaPHhwyCdG3nLLLXr44YcjfsBFsOBjjDGaMmWKRo0aVe5XyxyLttlZVyTffPONnn/+eb366qsxfUrpDz/8oMTExIj7ExIS9NNPP0Vd3zvvvBN2+4EDB7R69Wr95S9/0UsvvaQ//elPx7Rtdsdp5zyQpBdffFErVqwI+6FDS5cu1SWXXKKXX345qk/Ftms9Cm7b8uXLbWmbnXFWNlOBb4mzez0q2ZaKtC2a+mNh51wo2Y7MzExNmTIl7rbZWd/cuXPDPpbGs7bZ3bb7779fu3fvDruvdu3aev/997V27dqo6tq8eXPI1xWlpKSEfJXMGWecoZtuuinqttlZn91rm92x2jkX7F5DPvzww5D6/OcgwV+ZFS27zyvtbNtxw5k8GHaw65kst9ttxowZYzZv3hyyvSKvkKakpJT5vrnt27ebGjVqHPO2WZZlnnrqKVsv8WzevLlZsmSJMcZ3yZFlWea2226LuW3GGNOgQQPzxRdfGGN8lxZZlhXywURZWVmmTp06cdVdkVdK7G5bZcV58OBB8+yzz5rf//73xu12m7PPPjvmy4tat24dMq9Kevvtt6P+YI1oPPPMM+aMM86I6lg722Z3nHbOA2N8H5gyZcqUiPsnT55s+vXrF1Vddq1HldE2O+syxvcqU/Cz9bVq1Sr1Xrpolfdqa0U+GKai65HdbbOzPjvngp2veFdGfWWJZW0zxt622XkFWUpKitm4cWPE/VlZWSY5OTnqttlZn91rm92x2v24EMzusRtvfZVxXmlX244nJKTHsbS0tJAPdWjYsKHJzs4O/L5p0yaTmppabj0rV640I0eONKmpqeaMM84wTz/9tPnxxx8rlJCmpaWZDz74IOL+999/36SlpR3ztlXGtf4JCQnm+++/D/xeo0aNqD6QJ5wT5QTQ7jhXrFhhrrnmGlOrVi3z29/+1rjd7rg/Yfrmm282Xbp0CXuycvjwYdOlSxdzyy23xFV3OJs3bzZ169Y95m2zO04754ExvjWkrE/3/vzzz6NaQ/x12bEeVVbb7KrLGN/cGjhwoLn00kvNpZdeahISEky/fv0Cv/tv0SgvuXVyPbK7bXbWZ+dcsPMJhsqoryyxrG3G2Ns2O98L2bZt2zI/JPLNN980bdq0ibptdtZn99pmd6x2Py4Eq0oJaWV/NsmJkJByye5xLC8vL+T7xEpeXuf1elVQUFBuPRkZGcrIyNBTTz2lf/3rX5ozZ47Gjh0rr9erJUuWKD09PeS7nqLRq1cvPf300zrvvPPC7p8+fXqp7wQ8Fm2L5QvEo+X1ekMufXS73TF9uXhJJdtYGW2Ol51ts6OuRx99VHPmzNHPP/+sK664Qh9//LF+97vfKTExMabvbgw2ceJEzZ07V+3bt9fNN9+sDh06yLIsZWVlacaMGfJ4PJowYUJcdYfzyy+/KCUl5Zi3ze447Z4H+/bti/iF5ZKUlpam/fv3R1WXXetRZbTNzrqk0t/NG+n7b6NhSnynacnvM43m8aWy2N02O+uzcy6U1y6/aL5XtjLqK0ssa5vdbXvttddKvV3p9ddfV+vWrSVJr776qmbMmFHuW5ok31uH7r33Xg0aNKhUPL/88ovuv/9+DRo0qNx6KqM+u9c2u2O1+3GhKqpK52jHM8sYm994gWOmXbt2evjhhzVkyJCw+//f//t/Gj9+vL799tuY6960aZOef/75wJdK9+3bV/PmzYv6/uvWrVNGRoYGDx6su+66K/BFxhs3btSjjz6qBQsWaOXKlTF/wXJF2+ZyuZSbmxvyJeUV5XK5NGDAgMCD6LvvvqvzzjsvrgfR8uoqKCjQwoUL5fF4Ym5n7dq19cUXXwQekGNlZ9vsqishIUF///vfNWnSJLnd7sD2xMREffHFF+rcuXNcse7YsUM33nijFi1aFHhfk2VZ6t+/v2bOnKmWLVvGVW84t9xyi7Zu3arMzMxj3jY767JzHki+E5fc3Fw1atQo7P4ffvhBTZs2jWq82b0e2dk2O+uy2/Dhw6M67oUXXoi57oquR3a3zc767JwLVTnO8sS6ttnZtiZNmuiDDz7QKaecIklq1KiRPvvss8CatnnzZp1++unKy8srt64ffvhBXbt2VVJSkm6++ebA+wU3btyoZ555RsXFxVq3bl2ZTyxVVn12r212x2r340Kwiq4h4erbsGGDWrVqFdP9KuO8sqR423Y8ISE9jt12222BN+aHeyarR48euuCCC/TUU0/F/Tc8Ho/effddzZkzJ6aEVJLmz5+v6667Tnv37g3Z3qBBAz333HO66KKL4m5XRdtmJzsfRO2sa/r06SG///3vf9edd96phg0bhmy/9dZbo/qbVTHOhx56SC+++KKOHDmiK664QldddZW6dOlS4YTUb//+/fr2229ljFG7du3ietV17NixYbfn5eVpzZo12rp1q1asWKFu3bod87bZWZfdJ7olT2RKivXJGTvXIzvbZnecVZXd61FVdiyTPidV1tpmhxo1amj9+vWBBK2kjRs3qmvXrjpy5EhU9WVnZ+vGG2/UkiVLQp6869u3r2bOnBlzUmRnfXafa9nZtqp8TvPHP/4x5Hc7k+WKqsptqywkpMcxu5/Jqgy//PKLFi5cGDjZbd++vfr166eaNWs61qYTRTTPpFmWpW3bth2D1lSuZcuWac6cOXr77bfVpk0bff3111q2bJnOPvtsp5tW6pNT/VJTU9WxY0eNHj1aLVq0OMatqvoq46TervWoKj45U9WdSOvRiaIqr22VdQXZvn37Avdp27Zt3J/wand9lXGuZXesFWX3GlKV196q3LZKc+zerorKsG3bNtO/f//Adyb5v7uof//+jr4B+oMPPjCdOnWK+GXNnTt3NsuXL3egZajO8vLyzKxZs8zpp59uXC6XycjIME888YTTzYLDWI+AE8utt95qOnfuHPFD2zp37mxuvfVWB1pmL9Y2VBeu8lNWVGWtWrXSwoUL9dNPP+nTTz/Vp59+qp9++kkLFy607br6eEybNk3XX3+9UlNTS+2rU6eO/vrXv+rJJ590oGUnjqVLl6pz587Kz88vtS8vL0+nnHKKVqxY4UDLKk9qaqpGjRql1atX64svvtCZZ56phx9+2OlmwWGsR847EdcjOGf8+PHat2+fOnTooMcee0z//e9/NW/ePD366KPq0KGD9u/fr/HjxzvdzAo7kdY21pDqjUt2USlatGihhQsXqlOnTmH3b9y4Uf369dPOnTuPcctOHBdddJH69OkT8VMEp0+frg8//FDvvPPOMW6ZfUq+zyIcy7LUrFkz9e3bV3/4wx+OQatQ1bAeOe9EWI9Qtdj9vs+q6ERa21hDqjdeIUWl+OGHH0I+6rukhISEUl9TA3t98cUXuvDCCyPu79evn9auXXsMW2S/OnXqlHurVauWtmzZomHDhunee+91uslwAOuR806E9QhVS1W9gsxOJ9LaxhpSvfE9pKgUzZo105dffqm2bduG3b9hwwadfPLJx7hVJ5YT4YEqljf0L1iwQDfeeKMmTZpUiS1CVcR65LwTYT1C1VS/fn2dccYZTjejUpxIaxtrSPXGK6SoFP4vVw73keq//PKL7rvvPg0ePNiBlp04/A9UkVSnB6ponH322erRo4fTzYADWI+cx3oE2O9EWttYQ6o33kOKSvHDDz/otNNOk9vt1s0336wOHTrIsixlZWVpxowZ8ng8+vzzzx39Sprq7pZbbtFHH32kzz77LOz31J5xxhnq06dPqe/2Aqob1iPnsR4B9juR1jbWkOqNhBSVZseOHbrxxhu1aNGikA8U6N+/v2bOnKmWLVs628Bq7kR6oALKw3rkLNYjoHKcKGsba0j1RkKKSrd///7AlzW3a9dO9erVc7pJJ4wT5YEKiBbrkXNYj4DKcyKsbawh1RcJKXACOBEeqAAcH1iPAFQEa0j1Q0IKAAAAAHAEn7ILAAAAAHAECSkAAAAAwBEkpAAAAAAAR5CQAgAAAAAcQUIKAEAcdu3apREjRqhp06ZKSkpSixYtdNttt2nv3r1R17F9+3ZZlqX169dXXkMBAKjCSEgBAIjRtm3b1KNHD23evFlvvPGGvv32W82ePVsffPCBMjIytG/fPqebCADAcYGEFACAGN10001KSkrS4sWL1bt3bzVv3lwDBgzQ+++/r++//14TJkyQ5PvS9v/85z8h961bt65efPFFSVKrVq0kSd26dZNlWTr33HMDx82ZM0ennHKKkpOTdfLJJ+vmm28O7Nu5c6cuvvhi1apVS6mpqbrsssv0ww8/BPb/4x//UNeuXTVnzhw1b95ctWrV0o033iiPx6NHH31UTZo0UePGjTV58uSQtuXl5emGG25Q48aNlZqaqvPOO09ffPFFYP8XX3yhPn36qHbt2kpNTVX37t21Zs0aO/6lAIATFAkpAAAx2LdvnxYtWqTRo0erRo0aIfuaNGmiK6+8Um+++aai+Zrv1atXS5Lef/995eTkaO7cuZKkWbNm6aabbtINN9ygL7/8UvPmzVPbtm0lScYYXXLJJdq3b5+WLVumJUuWaOvWrRo2bFhI3Vu3btV7772nhQsX6o033tCcOXM0aNAgfffdd1q2bJkeeeQRTZw4UZ9++mmg3kGDBik3N1eZmZlau3atTjvtNJ1//vmBV3yvvPJK/eY3v9Fnn32mtWvX6u6771ZiYmLF/qEAgBNagtMNAADgeLJlyxYZY9SpU6ew+zt16qT9+/frp59+KreuRo0aSZIaNGigJk2aBLY/+OCD+tvf/qbbbrstsO3000+X5EteN2zYoOzsbKWnp0uSXnnlFZ1yyin67LPPAsd5vV7NmTNHtWvXVufOndWnTx9t2rRJmZmZcrlc6tChgx555BF99NFH+v3vf68PP/xQX375pX788UclJydLkh5//HH95z//0b///W/dcMMN2rlzp+6880517NhRktSuXbtY/30AAITgFVIAAGzkf2XUsqy47v/jjz9q9+7dOv/888Puz8rKUnp6eiAZlaTOnTurbt26ysrKCmxr2bKlateuHfg9LS1NnTt3lsvlCtn2448/SpLWrl2rn3/+WQ0aNFCtWrUCt+zsbG3dulWSNHbsWI0cOVIXXHCBHn744cB2AADiRUIKAEAM2rZtK8uy9M0334Tdv3HjRtWrV08NGzaUZVmlLt0tKioqs/6SlwGXZIwJm+yW3F7yUlrLssJu83q9knyvqJ588slav359yG3Tpk268847Jfnem/r1119r0KBBWrp0qTp37qx33nmnzPYCAFAWElIAAGLQoEED9e3bVzNnztQvv/wSsi83N1evvfaahg0bJsuy1KhRI+Xk5AT2b9myRYcPHw78npSUJEnyeDyBbbVr11bLli31wQcfhP37nTt31s6dO7Vr167Atm+++UZ5eXkRLyOOxmmnnabc3FwlJCSobdu2IbeGDRsGjmvfvr3GjBmjxYsX649//KNeeOGFuP8mAAAkpAAAxOiZZ55RQUGB+vfvr+XLl2vXrl1auHCh+vbtq2bNmgU+vfa8887TM888o88//1xr1qzRqFGjQl6lbNy4sWrUqKGFCxfqhx9+UF5eniTfK5FPPPGEpk+fri1btujzzz/X008/LUm64IILdOqpp+rKK6/U559/rtWrV+vqq69W79691aNHj7hjuuCCC5SRkaFLLrlEixYt0vbt27Vy5UpNnDhRa9as0S+//KKbb75ZH330kXbs2KFPPvlEn332WYWSYAAASEgBAIhRu3bttGbNGrVp00bDhg1TmzZtdMMNN6hPnz5atWqV6tevL0l64oknlJ6erl69eunPf/6z7rjjDtWsWTNQT0JCgqZPn67/+7//U9OmTXXxxRdLkq655hpNmzZNM2fO1CmnnKLBgwdry5Ytko5+lUy9evXUq1cvXXDBBWrdurXefPPNCsVkWZYyMzPVq1cvXXfddWrfvr0uv/xybd++XWlpaXK73dq7d6+uvvpqtW/fXpdddpkGDBig+++/v0J/FwBwYrNMNJ9LDwAAAACAzXiFFAAAAADgCBJSAAAAAIAjSEgBAAAAAI4gIQUAAAAAOIKEFAAAAADgCBJSAAAAAIAjSEgBAAAAAI4gIQUAAAAAOIKEFAAAAADgCBJSAAAAAIAjSEgBAAAAAI4gIQUAAAAAOOL/AzTjrOAywAqIAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1100x700 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "monopoly_barchart(P)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "\n",
    "\n",
    "\n",
    "With 40 squares, the average is 2.5% (shown with the dotted red line), but a few squares stand out as deviating from the average:  `JAIL` is  over 6% (due to the \"Go to Jail\" space, the two \"Go to Jail\" cards, and the three-doubles rule). `E3` (Illinois Avenue) is the next most probable square, helped by a community chest card and by the fact that it is 14 spaces from jail, meaning it is the most probable place to land after two rolls out from jail. `G2J` (\"Go to Jail\") is 0%, because you can't end a roll there. <img src=\"https://upload.wikimedia.org/wikipedia/commons/thumb/f/f1/Dice_simple_flat.svg/239px-Dice_simple_flat.svg.png\" style=\"float:right\">After that square, the three chance squares (`CH1`, `CH2`, and `CH3`) are the least probable, because 10 of the 16 chance cards send the player away from the square. The first three railroads are near 3%, helped by the chance cards, but `R4` (Short Line) is only at about 2.5%, because no chance card sends you there.\n",
    "\n",
    "If you'd like to see another simulation, check out my [Dice Baseball](Dice%20Baseball.ipynb) notebook, or the example below:\n",
    "\n",
    "# Simulating another Game: Roll-1\n",
    "\n",
    "\n",
    "To better understand how simulations work here is another problem: a dice game (due to my friend [Cam Hewett](https://www.mathgenealogy.org/id.php?id=274478)) that I call *Roll-1*. Two players take turns rolling a single die, and the first to roll a 1 wins. The complication is that we start with a *d*-sided die, and whenever a player rolls any number other than 1, the next player's roll is with a die with that number of sides.  The function `roll_1(d)` simulates a single game (starting with a *d*-sided die) and returns the number of rolls in the game:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def roll_1(d=6) -> int:\n",
    "    \"\"\"Simulate one game of Roll-1, starting with a d-sided die; return number of rolls in the game.\"\"\"\n",
    "    for turn in count_from(1):\n",
    "        d = roll(d)\n",
    "        if d == 1:\n",
    "            return turn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function `monte_carlo` applies the Monte Carlo method, calling a simulation function *N* times and collecting a probability distribution of the results. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def monte_carlo(simulation: Callable, *args, N=1_000_000) -> ProbDist:\n",
    "    \"\"\"Monte Carlo repetition of simulation(*args) for N trials, as a ProbDist of outcomes.\"\"\"\n",
    "    return ProbDist(simulation(*args) for _ in range(N))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is the full probability distribution of rolls (starting with a 6-sided die):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGxCAYAAACKvAkXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA8CklEQVR4nO3de3wU9b3/8ffs5oaVRK4JYBIQ5RJRgeAlUMALBClYsFpjtYAKeihohbT9VYQWxQvqqRpRQW3VVI8iD49SPbpKYqlECRYJSfESbVqCiZJIQ0oS0Fx2Z35/bDNmyYUkbLLAvJ6Pxz7yyXdnJt/PzATfzk52DcuyLAEAADiIK9QTAAAA6G4EIAAA4DgEIAAA4DgEIAAA4DgEIAAA4DgEIAAA4DgEIAAA4DgEIAAA4DhhoZ7Ascg0Te3du1c9e/aUYRihng4AAGgHy7JUU1OjgQMHyuVq+xoPAagFe/fuVXx8fKinAQAAOqG0tFSnnnpqm8uEPACtXbtW//3f/62ysjKdeeaZysjI0MSJE1tc9tVXX9W6detUUFCguro6nXnmmbrjjjs0bdo0e5nMzExdf/31zdb99ttvFRUV1a459ezZU5J/B0ZHR3eiKwAA0N2qq6sVHx9v/3e8LSENQBs2bNCSJUu0du1aTZgwQU8++aSmT5+uTz/9VAkJCc2Wz8nJ0dSpU3XvvffqlFNO0bPPPqvLLrtMf/3rXzVmzBh7uejoaH3++ecB67Y3/EiyX/aKjo4mAAEAcJxpz+0rRig/DPX888/X2LFjtW7dOnts5MiRmj17tlavXt2ubZx55plKS0vTb3/7W0n+K0BLlizRgQMHOj2v6upqxcTEqKqqigAEAMBxoiP//Q7ZX4HV19crLy9PqampAeOpqanKzc1t1zZM01RNTY169+4dMH7w4EElJibq1FNP1cyZM5Wfn9/mdurq6lRdXR3wkCSfz2d/ban2er0BtWmabdYNDQ0BdWP2bKwty2pWSwqoTdMMqL1eb5u1z+cLqOmJnuiJnuiJnk7kntorZAGooqJCPp9PsbGxAeOxsbEqLy9v1zYefPBBHTp0SFdddZU9NmLECGVmZur111/X+vXrFRUVpQkTJqioqKjV7axevVoxMTH2o/EG6I8//liSVFhYqMLCQknSrl277G3l5+eruLhYkrR9+3aVlpZKknJzc1VWVibJ/7JdRUWFJGnz5s32lamsrCzV1NRIkjwej2pra+X1euXxeOT1elVbWyuPxyNJqqmpUVZWliTpwIED2rx5s70Pc3JyJEllZWV2cCwtLdX27dslScXFxXYALCoq0q5du+iJnuiJnuiJnk7InrZu3ap2s0Lkq6++siRZubm5AeN33323NXz48COu/+KLL1onnXSSlZ2d3eZyPp/POuecc6xbbrml1WVqa2utqqoq+1FaWmpJsiorKy3Lsiyv12t5vd5mdUNDQ0Dt8/narOvr6wNq0zQDatM0m9WWZQXUPp8voG5oaGiz9nq9AXVLfdATPdETPdETPZ0IPe3fv9+SZFVVVVlHErJ7gOrr63XSSSfp5Zdf1uWXX26P33rrrSooKNCWLVtaXXfDhg26/vrr9fLLL2vGjBlH/Fk33nijvvzyS7311lvtmhv3AAEAcPw5Lu4BioiIUHJysrKzswPGs7OzNX78+FbXW79+va677jq9+OKL7Qo/lmWpoKBAAwYMOOo5AwCAE0NI/ww+PT1dc+bM0bhx45SSkqKnnnpKJSUlWrhwoSRp2bJl+uqrr/Tcc89J8oefuXPn6pFHHtEFF1xg3yvUo0cPxcTESJLuvPNOXXDBBTrjjDNUXV2tNWvWqKCgQI8//nhomgQAAMeckAagtLQ07d+/X6tWrVJZWZlGjRolj8ejxMRESf4bp0pKSuzln3zySXm9Xi1evFiLFy+2x+fNm6fMzExJ/puwbrrpJpWXlysmJkZjxoxRTk6OzjvvvG7tDQAAHLtC+j5AxyruAQIA4PhzXNwDBAAAECoEIAAA4DgEIAAA4DgEIAAA4DgEIAAA4DgEIAAA4DghfR8gHL3Bt73Z6XX33Hfkd9IGAOBExBUgAADgOAQgAADgOAQgAADgOAQgAADgOAQgAADgOAQgAADgOAQgAADgOAQgAADgOAQgAADgOAQgAADgOAQgAADgOAQgAADgOAQgAADgOAQgAADgOAQgAADgOAQgAADgOAQgAADgOAQgAADgOAQgAADgOAQgAADgOAQgAADgOAQgAADgOAQgAADgOAQgAADgOAQgAADgOAQgAADgOAQgAADgOAQgAADgOAQgAADgOAQgAADgOAQgAADgOAQgAADgOAQgAADgOAQgAADgOAQgAADgOAQgAADgOGGhngCOHYNve7PT6+65b0YQZwIAQNfiChAAAHAcAhAAAHAcAhAAAHAcAhAAAHAcAhAAAHAcAhAAAHAcAhAAAHAcAhAAAHAcAhAAAHAcAhAAAHAcAhAAAHAcAhAAAHAcAhAAAHAcAhAAAHAcAhAAAHAcAhAAAHAcAhAAAHAcAhAAAHCckAegtWvXasiQIYqKilJycrLee++9Vpd99dVXNXXqVPXr10/R0dFKSUnRpk2bmi33yiuvKCkpSZGRkUpKStLGjRu7sgUAAHCcCWkA2rBhg5YsWaLly5crPz9fEydO1PTp01VSUtLi8jk5OZo6dao8Ho/y8vJ00UUX6bLLLlN+fr69zLZt25SWlqY5c+bob3/7m+bMmaOrrrpKf/3rX7urLQAAcIwzLMuyQvXDzz//fI0dO1br1q2zx0aOHKnZs2dr9erV7drGmWeeqbS0NP32t7+VJKWlpam6ulpvvfWWvcyll16qXr16af369e3aZnV1tWJiYlRVVaXo6OgOdNT9Bt/2ZqfX3XPfjC7bFgAA3a0j//0O2RWg+vp65eXlKTU1NWA8NTVVubm57dqGaZqqqalR79697bFt27Y12+a0adPa3GZdXZ2qq6sDHpLk8/nsry3VXq83oDZNs826oaEhoG7Mno21ZVnNakkBtWmaAXWEy78Nl2HZtfuwOryFOsywWuwp3GXJbfiXiTisdjWt5a8jm9TB6snr9bZZ+3y+gPp4OE70RE/0RE/01D09tVfIAlBFRYV8Pp9iY2MDxmNjY1VeXt6ubTz44IM6dOiQrrrqKnusvLy8w9tcvXq1YmJi7Ed8fLwk6eOPP5YkFRYWqrCwUJK0a9cuFRUVSZLy8/NVXFwsSdq+fbtKS0slSbm5uSorK5Pkf9muoqJCkrR582YdOHBAkpSVlaWamhpJksfjUW1trbxerzwej7xer2pra+XxeCRJNTU1ysrKkiQdOHBAmzdvtvdh+ln+k+HsXpYWJfnrc/tZumG4/2SYGGfp2tP99dRBpq4Y7K9nJpgt9nTt6aYmxvlP/BuGmzq3n79elOTT2b38dfpZPg2L8de3j/Yp4WQFtaecnBxJUllZmR1cS0tLtX37dklScXGx/bJnUVGRdu3adcwfJ3qiJ3qiJ3rq+p62bt2q9grZS2B79+7VoEGDlJubq5SUFHv8nnvu0fPPP6/PPvuszfXXr1+vBQsW6LXXXtOUKVPs8YiICP3xj3/UT37yE3vshRde0Pz581VbW9viturq6lRXV2d/X11drfj4eFVWVqpXr1524nS73QG11+uVYRh27XK55HK5Wq0bGhrkdrvtOiwsTIZh2LXkT69N6/DwcFmWZdemacrn89n1iBUe1ZuGXIalMEOqNw25DUvuJrXLkBoOq8MMS5/fPT2gp6HL31a4y5JpST7LUITLkq9J7bUks7E2JVOGIl2WGkxp930zg9aTaZoKCwtrtfb5fLIsy65bOjbH2nGiJ3qiJ3qip67vqbKyUn369GnXS2BhbT7bhfr27Su3293sysy+ffuaXcE53IYNGzR//ny9/PLLAeFHkuLi4jq8zcjISEVGRjYbd7vdAV8PrxtPhPbW4eHhnaoNw7DrxoPcWNebhiR/MKn/T5T1WYZ8R6i9ltFifw3/2Z4ke9tt1XVN6mD1dKS6teNxLB8neqIneqIneureno4kZC+BRUREKDk5WdnZ2QHj2dnZGj9+fKvrrV+/Xtddd51efPFFzZjR/MbblJSUZtvMyspqc5sAAMBZQnYFSJLS09M1Z84cjRs3TikpKXrqqadUUlKihQsXSpKWLVumr776Ss8995wkf/iZO3euHnnkEV1wwQX2lZ4ePXooJiZGknTrrbdq0qRJuv/++zVr1iy99tpreuedd/T++++HpkkAAHDMCen7AKWlpSkjI0OrVq3S6NGjlZOTI4/Ho8TEREn+G6eavifQk08+Ka/Xq8WLF2vAgAH249Zbb7WXGT9+vF566SU9++yzOvvss5WZmakNGzbo/PPP7/b+AADAsSmk7wN0rOJ9gI5+WwAAdLfj4n2AAAAAQoUABAAAHIcABAAAHIcABAAAHIcABAAAHIcABAAAHIcABAAAHIcABAAAHIcABAAAHIcABAAAHIcABAAAHIcABAAAHIcABAAAHIcABAAAHIcABAAAHIcABAAAHIcABAAAHIcABAAAHIcABAAAHIcABAAAHIcABAAAHIcABAAAHIcABAAAHIcABAAAHIcABAAAHIcABAAAHIcABAAAHIcABAAAHIcABAAAHIcABAAAHIcABAAAHIcABAAAHIcABAAAHIcABAAAHIcABAAAHIcABAAAHIcABAAAHIcABAAAHIcABAAAHIcABAAAHIcABAAAHIcABAAAHCcs1BNwosG3vdnpdffcNyOIMwEAwJm4AgQAAByHAAQAAByHAAQAAByHAAQAAByHAAQAAByHAAQAAByHAAQAAByHAAQAAByHAAQAAByHAAQAAByHAAQAAByHAAQAAByHAAQAAByHAAQAAByHAAQAAByHAAQAABwnLNQTwIlp8G1vdnrdPffNCOJMAABojitAAADAcQhAAADAcUIegNauXashQ4YoKipKycnJeu+991pdtqysTNdcc42GDx8ul8ulJUuWNFsmMzNThmE0e9TW1nZhFwAA4HgS0gC0YcMGLVmyRMuXL1d+fr4mTpyo6dOnq6SkpMXl6+rq1K9fPy1fvlznnHNOq9uNjo5WWVlZwCMqKqqr2gAAAMeZkAaghx56SPPnz9eCBQs0cuRIZWRkKD4+XuvWrWtx+cGDB+uRRx7R3LlzFRMT0+p2DcNQXFxcwAMAAKBRyAJQfX298vLylJqaGjCempqq3Nzco9r2wYMHlZiYqFNPPVUzZ85Ufn5+m8vX1dWpuro64CFJPp/P/tpS7fV6A2rTNNusGxoa7DrSbcmQdVhtKdLt/2rYtQJq0zTV0NBg1xEu/7jLsOzafVgd3kIdZlgt9hTusuQ2/MtEHFa7mtaNc29SNzQ0yLI63pNLliJd39Ver9fur6Xa5/MF1F15nJr21FhbltWslhRQH36c6Ime6Ime6Kl7emqvTgWgzMxMffPNN51Z1VZRUSGfz6fY2NiA8djYWJWXl3d6uyNGjFBmZqZef/11rV+/XlFRUZowYYKKiopaXWf16tWKiYmxH/Hx8ZKkjz/+WJJUWFiowsJCSdKuXbvsbeXn56u4uFiStH37dpWWlkqScnNzVVZWJknKyclRRUWFJGnz5s06cOCAJGnVWJ9ie/h//gPn+RQdIUW6/XWkW4qO8NeSFNvDv7wkHThwQJs3b5bk34fpZ/nHz+5laVGSvz63n6UbhvtPholxlq493V9PHWTqisH+emaC2WJP155uamKc/8S/Ybipc/v560VJPp3dy1+nn+XTsBh/fftonxJO9veRlZWlmpqaDveUcLJ/O5I0LMZSTk6OJP89X41huLS0VNu3b5ckFRcX26G2qKhIu3bt6rLj1LQnj8ej2tpaeb1eeTweeb1e1dbWyuPxSJJqamqUlZXV4nGiJ3qiJ3qip67vaevWrWovw2qMeR0wYMAAHTp0SD/+8Y81f/58jR8/vqOb0N69ezVo0CDl5uYqJSXFHr/nnnv0/PPP67PPPmtz/QsvvFCjR49WRkZGm8uZpqmxY8dq0qRJWrNmTYvL1NXVqa6uzv6+urpa8fHxqqysVK9evezE6Xa7A2qv1yvDMOza5XLJ5XK1Wjc0NMjtduu0299SpNtSvc9/XeS72h8Y6nySISnCLdX5DBmy7Hr3vdPl8/kUHh4u0zQ1YoVH9aYhl2EpzJDqTUNuw5K7Se0ypIbD6jDD0ud3Tw/oaejytxXusmRaks8yFOGy5GtSey3JbKxNyZShSJelBlPafd9MNTQ0KCwsTIZhaPjyN9rdk0uWwl1Snemv/373pQoLC5NpmjJNs1nt8/lkWZZdt3RsgnGcGuvGnhpryf9/GU3r8PBwWZZl16ZpBhynlvqgJ3qiJ3qip+D2VFlZqT59+qiqqkrR0dFqS6feCPHLL7/Um2++qczMTF100UUaMmSIrr/+es2bN6/d99v07dtXbre72dWeffv2NbsqdDRcLpfOPffcNq8ARUZGKjIystm42+0O+Hp43XgitLcODw+36zqf0Urt/2oF1IZdNx7kxrre9K9rWobq/xNlfZYh3xFqr2W02F+D+d1c6ttR1zWpW++v7Z5MGaozv6sb99nhvTbWrR2PrjhOHakNw7Dr1uZOT/RET/RET93T05F06iUwt9utH/7wh3r11VdVWlqqm266SS+88IISEhL0wx/+UK+99pr9elxrIiIilJycrOzs7IDx7OzsTl1Rao1lWSooKNCAAQOCtk0AAHB8O+qboPv3768JEyYoJSVFLpdLH330ka677joNHTpU7777bpvrpqen6w9/+IOeeeYZFRYWaunSpSopKdHChQslScuWLdPcuXMD1ikoKFBBQYEOHjyof/3rXyooKNCnn35qP3/nnXdq06ZN2r17twoKCjR//nwVFBTY2wQAAOj0Z4F9/fXXev755/Xss89q9+7dmj17tt544w1NmTJF3377rVasWKF58+bpiy++aHUbaWlp2r9/v1atWqWysjKNGjVKHo9HiYmJkvw3Th3+nkBjxoyx67y8PL344otKTEzUnj17JPlvwrrppptUXl6umJgYjRkzRjk5OTrvvPM62yoAADjBdCoAXXbZZdq0aZOGDRumG2+8UXPnzlXv3r3t53v06KFf/OIXevjhh4+4rUWLFmnRokUtPpeZmdls7Ej3bD/88MPt+rkAAMC5OhWA+vfvry1btgT89dbhBgwYYP+5GgAAwLGkU/cATZ48WWPHjm02Xl9fr+eee06S/87wxpeyAAAAjiWdCkDXX3+9qqqqmo3X1NTo+uuvP+pJAQAAdKVOBSDLsmQYRrPxL7/8ss3P6AIAADgWdOgeoDFjxsgwDBmGoUsuuSTgDYd8Pp+Ki4t16aWXBn2SAAAAwdShADR79mxJ/vfimTZtmk4++WT7uYiICA0ePFhXXHFFUCcIAAAQbB0KQCtXrpQkDR48WGlpaYqKiuqSSQEAAHSlTv0Z/Lx584I9DwAAgG7T7gDUu3dv/f3vf1ffvn3Vq1evFm+CblRZWRmUyQEAAHSFdgeghx9+WD179rTrtgIQAADAsazdAajpy17XXXddV8wFAACgW7Q7AFVXV7d7o9HR0Z2aDAAAQHdodwA65ZRTjviyV+MbJPp8vqOeGAAAQFdpdwD6y1/+0pXzAAAA6DbtDkCTJ0/uynkAAAB0m3YHoF27dmnUqFFyuVzatWtXm8ueffbZRz0xAACArtLuADR69GiVl5erf//+Gj16tAzDkGVZzZbjHiAAAHCsa3cAKi4uVr9+/ewaAADgeNXuAJSYmNhiDQAAcLzp1GeBSdLnn3+uRx99VIWFhTIMQyNGjNAtt9yi4cOHB3N+AAAAQefqzEr/+7//q1GjRikvL0/nnHOOzj77bO3cuVOjRo3Syy+/HOw5AgAABFWnrgD9v//3/7Rs2TKtWrUqYHzlypX69a9/rR//+MdBmRwAAEBX6NQVoPLycs2dO7fZ+E9/+lOVl5cf9aQAAAC6UqcC0IUXXqj33nuv2fj777+viRMnHvWkAAAAulK7XwJ7/fXX7fqHP/yhfv3rXysvL08XXHCBJOmDDz7Qyy+/rDvvvDP4swQAAAiidgeg2bNnNxtbu3at1q5dGzC2ePFiLVy48KgnBgAA0FXaHYBM0+zKeQAAAHSbTt0DBAAAcDzr9BshHjp0SFu2bFFJSYnq6+sDnvv5z39+1BMDAADoKp0KQPn5+frBD36gb775RocOHVLv3r1VUVGhk046Sf379ycAAQCAY1qnXgJbunSpLrvsMlVWVqpHjx764IMP9MUXXyg5OVm/+93vgj1HAACAoOpUACooKNAvfvELud1uud1u1dXVKT4+Xg888IBuv/32YM8RAAAgqDoVgMLDw2UYhiQpNjZWJSUlkqSYmBi7BgAAOFZ16h6gMWPGaMeOHRo2bJguuugi/fa3v1VFRYWef/55nXXWWcGeIwAAQFB16grQvffeqwEDBkiS7rrrLvXp00c/+9nPtG/fPj311FNBnSAAAECwdeoK0Lhx4+y6X79+8ng8QZsQAABAV+v0+wBJ0r59+/T555/LMAwNHz5c/fr1C9a8AAAAukynXgKrrq7WnDlzNGjQIE2ePFmTJk3SwIED9dOf/lRVVVXBniMAAEBQdSoALViwQH/961/1xhtv6MCBA6qqqtIbb7yhHTt26MYbbwz2HAEAAIKqUy+Bvfnmm9q0aZO+//3v22PTpk3T73//e1166aVBmxwAAEBX6NQVoD59+igmJqbZeExMjHr16nXUkwIAAOhKnQpAK1asUHp6usrKyuyx8vJy/epXv9JvfvOboE0OAACgK7T7JbAxY8bY7/4sSUVFRUpMTFRCQoIkqaSkRJGRkfrXv/6l//qv/wr+TAEAAIKk3QFo9uzZXTgNAACA7tPuALRy5cqunAcAAEC3Oao3QszLy1NhYaEMw1BSUpLGjBkTrHkBAAB0mU4FoH379unqq6/Wu+++q1NOOUWWZamqqkoXXXSRXnrpJd4RGgAAHNM69Vdgt9xyi6qrq/XJJ5+osrJS//73v/Xxxx+rurpaP//5z4M9RwAAgKDq1BWgt99+W++8845GjhxpjyUlJenxxx9Xampq0CYHAADQFTp1Bcg0TYWHhzcbDw8Pl2maRz0pAACArtSpAHTxxRfr1ltv1d69e+2xr776SkuXLtUll1wStMkBAAB0hU4FoMcee0w1NTUaPHiwhg4dqtNPP11DhgxRTU2NHn300WDPEQAAIKg6dQ9QfHy8du7cqezsbH322WeyLEtJSUmaMmVKsOcHAAAQdB0OQF6vV1FRUSooKNDUqVM1derUrpgXAABAl+nwS2BhYWFKTEyUz+frivkAAAB0uU5/GvyyZctUWVkZ7PkAAAB0uU7dA7RmzRr94x//0MCBA5WYmKjvfe97Ac/v3LkzKJMDAADoCp0KQLNnz5ZhGLIsK9jzAQAA6HIdCkDffPONfvWrX+lPf/qTGhoadMkll+jRRx9V3759u2p+AAAAQdehe4BWrlypzMxMzZgxQz/5yU/0zjvv6Gc/+1lXzQ0AAKBLdOgK0Kuvvqqnn35aV199tSTp2muv1YQJE+Tz+eR2u7tkggAAAMHWoStApaWlmjhxov39eeedp7CwsICPxAAAADjWdSgA+Xw+RUREBIyFhYXJ6/V2egJr167VkCFDFBUVpeTkZL333nutLltWVqZrrrlGw4cPl8vl0pIlS1pc7pVXXlFSUpIiIyOVlJSkjRs3dnp+AADgxNOhl8Asy9J1112nyMhIe6y2tlYLFy4M+FP4V199tV3b27Bhg5YsWaK1a9dqwoQJevLJJzV9+nR9+umnSkhIaLZ8XV2d+vXrp+XLl+vhhx9ucZvbtm1TWlqa7rrrLl1++eXauHGjrrrqKr3//vs6//zzO9IuAAA4QXXoCtC8efPUv39/xcTE2I+f/vSnGjhwYMBYez300EOaP3++FixYoJEjRyojI0Px8fFat25di8sPHjxYjzzyiObOndvqz8nIyNDUqVO1bNkyjRgxQsuWLdMll1yijIyMjrQKAABOYB26AvTss88G7QfX19crLy9Pt912W8B4amqqcnNzO73dbdu2aenSpQFj06ZNazMA1dXVqa6uzv6+urpakuyP+2j86na7A2qv1yvDMOza5XLJ5XK1Wjc0NNg3i0e6LdX7JEtGk1qKdEt1PsmQFOGW6nyGDFl2bZqmfD6fwsPDZZqmIlyW6k1DLsNSmCHVm4bchiV3k9plSA2H1WGGZd+83vRjTcJdlkxL8lmGIlyWfE1qryWZjbUpmTIU6bLUYPrXbWhoUFhYmAyjYz25ZCncJdWZ/trr9SosLEymaco0zWa1z+eTZVl23dKxCcZxaqwbe2qsJdlzbKzDw8NlWZZdH36cWuqDnuiJnuiJnoLfU3t16qMwgqGiokI+n0+xsbEB47GxsSovL+/0dsvLyzu8zdWrVwdcwYqPj5ckffzxx5KkwsJCFRYWSpJ27dqloqIiSVJ+fr6Ki4slSdu3b1dpaakkKTc3V2VlZZKknJwcVVRUSJI2b96sAwcOSJJWjfUptof/5z9wnk/REf6g8MB5PkW6pegIfy1JsT38y0vSgQMHtHnzZkn+fZh+ln/87F6WFiX563P7WbphuD+VTIyzdO3p/nrqIFNXDPbXMxPMFnu69nRTE+P8b3B5w3BT5/bz14uSfDq7l79OP8unYTH++vbRPiWc7O8jKytLNTU1He4p4WT/diRpWIylnJwcSf57vhrDcGlpqbZv3y5JKi4uVn5+viSpqKhIu3bt6rLj1LQnj8ej2tpaeb1eeTweeb1e1dbWyuPxSJJqamqUlZXV4nGiJ3qiJ3qip67vaevWrWovwwrR2znv3btXgwYNUm5urlJSUuzxe+65R88//7w+++yzNte/8MILNXr06GZXdiIiIvTHP/5RP/nJT+yxF154QfPnz1dtbW2L22rpClB8fLwqKyvVq1evoCfs025/q9NXgHbfOz0gYY9Y4en0FaDP754e0NPQ5W93+grQ7vtmBvxfw/Dlb3T6CtDf776U/xOiJ3qiJ3qipw73VFlZqT59+qiqqkrR0dFqS6c+CiMY+vbtK7fb3ezKzL59+5pdwemIuLi4Dm8zMjIy4MbuRo0vVzV9j6OmdeOJ0N46PDzcrut8Riu1/6sVUBt23XiQG+t607+uaRmq/0+U9VmGfEeovZbRYn8N5ndzqW9HXdekbr2/tnsyZajO/K5u3GeH99pYt3Y8uuI4daQ2DMOuW5s7PdETPdETPXVPT0cSsgAUERGh5ORkZWdn6/LLL7fHs7OzNWvWrE5vNyUlRdnZ2QH3AWVlZWn8+PFHNV+EzuDb3uz0unvumxHEmQAAThQhC0CSlJ6erjlz5mjcuHFKSUnRU089pZKSEi1cuFCStGzZMn311Vd67rnn7HUKCgokSQcPHtS//vUvFRQUKCIiQklJSZKkW2+9VZMmTdL999+vWbNm6bXXXtM777yj999/v9v7AwAAx6aQBqC0tDTt379fq1atUllZmUaNGiWPx6PExERJ/hunSkpKAtYZM2aMXefl5enFF19UYmKi9uzZI0kaP368XnrpJa1YsUK/+c1vNHToUG3YsIH3AAIAALaQBiBJWrRokRYtWtTic5mZmc3G2nPP9pVXXqkrr7zyaKcGAABOUCH7M3gAAIBQIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHIQABAADHCQv1BIDuNPi2Nzu97p77ZgRxJgCAUOIKEAAAcBwCEAAAcBwCEAAAcBwCEAAAcBwCEAAAcBwCEAAAcJyQB6C1a9dqyJAhioqKUnJyst577702l9+yZYuSk5MVFRWl0047TU888UTA85mZmTIMo9mjtra2K9sAAADHkZAGoA0bNmjJkiVavny58vPzNXHiRE2fPl0lJSUtLl9cXKwf/OAHmjhxovLz83X77bfr5z//uV555ZWA5aKjo1VWVhbwiIqK6o6WAADAcSCkb4T40EMPaf78+VqwYIEkKSMjQ5s2bdK6deu0evXqZss/8cQTSkhIUEZGhiRp5MiR2rFjh373u9/piiuusJczDENxcXHd0gMAADj+hOwKUH19vfLy8pSamhownpqaqtzc3BbX2bZtW7Plp02bph07dqihocEeO3jwoBITE3Xqqadq5syZys/Pb3MudXV1qq6uDnhIks/ns7+2VHu93oDaNM0264aGBruOdFsyZB1WW4p0+78adq2A2jRNu1fTNBXh8o+7DMuu3YfV4S3UYYbVYk/hLktuw79MxGG1q2ndOPcmdUNDgyyr4z25ZCnS9V3t9Xrt/hrr9vYU9p85Nu2j6XHqaE9Nj1Njf5ZlNaslBdSHH6eWempa+3y+gLorz72mx4me6Ime6OlE66m9QhaAKioq5PP5FBsbGzAeGxur8vLyFtcpLy9vcXmv16uKigpJ0ogRI5SZmanXX39d69evV1RUlCZMmKCioqJW57J69WrFxMTYj/j4eEnSxx9/LEkqLCxUYWGhJGnXrl32tvLz81VcXCxJ2r59u0pLSyVJubm5KisrkyTl5OTYc9u8ebMOHDggSVo11qfYHv6f/8B5PkVHSJFufx3plqIj/LUkxfbwLy9JBw4c0ObNm+19mH6Wf/zsXpYWJfnrc/tZumG4/2SYGGfp2tP99dRBpq4Y7K9nJpgt9nTt6aYmxvlP/BuGmzq3n79elOTT2b38dfpZPg2L8de3j/Yp4WR/H1lZWaqpqelwTwkn+7cjScNiLOXk5EiSysrK7DDc3p5mJvjrKwabmjrIbHacOtpT0+Pk8XhUW1srr9crj8cjr9er2tpaeTweSVJNTY2ysrJaPE4t9VRaWqrt27dL8r+82xjUi4qKtGvXLkldc+41PU70RE/0RE8nUk9bt25VexlWY8zrZnv37tWgQYOUm5urlJQUe/yee+7R888/r88++6zZOsOGDdP111+vZcuW2WNbt27V97//fZWVlbX4spdpmho7dqwmTZqkNWvWtDiXuro61dXV2d9XV1crPj5elZWV6tWrl5043W53QO31emUYhl27XC65XK5W64aGBrndbp12+1uKdFuq9/mvi3xX+wNDnU8yJEW4pTqfIUOWXe++d7p8Pp/Cw8NlmqZGrPCo3jTkMiyFGVK9achtWHI3qV2G1HBYHWZY+vzu6QE9DV3+tsJdlkxL8lmGIlyWfE1qryWZjbUpmTIU6bLUYEq775uphoYGhYWFyTAMDV/+Rrt7cslSuEuqM/313+++VGFhYTJNU6Zp6vQVm9rdkyR5LcPu45+rZwYcp2G3v9HunsyAY2Oo6K5UhYX5XzX2er0BdXh4uCzLsmvTNAOOk2maAT0dXvt8PlmWZdctnW/BOPca68bj1FjTEz3REz2dCD1VVlaqT58+qqqqUnR0tNoSsnuA+vbtK7fb3exqz759+5pd5WkUFxfX4vJhYWHq06dPi+u4XC6de+65bV4BioyMVGRkZLNxt9sd8PXwuvFEaG8dHh5u13U+o5Xa/9UKqA27bjzIjXW96V/XtAzV/yfK+ixDviPUXstosb8G87u51LejrmtSt95f2z2ZMlRnflc37rOmvba3p5b6aHoMOtpT0z6a9tdSbRiGXR9+nI5Ut3aOdcW515GanuiJnujpeO3pSEL2ElhERISSk5OVnZ0dMJ6dna3x48e3uE5KSkqz5bOysjRu3LiAA9OUZVkqKCjQgAEDgjNxAABw3Avpn8Gnp6frD3/4g5555hkVFhZq6dKlKikp0cKFCyVJy5Yt09y5c+3lFy5cqC+++ELp6ekqLCzUM888o6efflq//OUv7WXuvPNObdq0Sbt371ZBQYHmz5+vgoICe5sAAAAh/TP4tLQ07d+/X6tWrVJZWZlGjRolj8ejxMRESf4bp5q+J9CQIUPk8Xi0dOlSPf744xo4cKDWrFkT8CfwBw4c0E033aTy8nLFxMRozJgxysnJ0Xnnndft/QEAgGNTSAOQJC1atEiLFi1q8bnMzMxmY5MnT9bOnTtb3d7DDz+shx9+OFjTAwAAJ6CQfxQGAABAdyMAAQAAxyEAAQAAxyEAAQAAxyEAAQAAxyEAAQAAxyEAAQAAxyEAAQAAxyEAAQAAxyEAAQAAxyEAAQAAxyEAAQAAxyEAAQAAxyEAAQAAxyEAAQAAxyEAAQAAxyEAAQAAxyEAAQAAxyEAAQAAxyEAAQAAxwkL9QSA49Hg2948qvX33DcjSDMBAHQGV4AAAIDjEIAAAIDjEIAAAIDjEIAAAIDjEIAAAIDjEIAAAIDjEIAAAIDjEIAAAIDjEIAAAIDjEIAAAIDjEIAAAIDjEIAAAIDjEIAAAIDjEIAAAIDjEIAAAIDjEIAAAIDjEIAAAIDjEIAAAIDjEIAAAIDjEIAAAIDjhIV6AoDTDb7tzaNaf899M4I0EwBwDq4AAQAAxyEAAQAAxyEAAQAAxyEAAQAAxyEAAQAAxyEAAQAAxyEAAQAAxyEAAQAAxyEAAQAAxyEAAQAAxyEAAQAAx+GzwIATCJ8rBgDtwxUgAADgOAQgAADgOAQgAADgOAQgAADgOAQgAADgOPwVGIAWHc1flPHXZACOdVwBAgAAjhPyALR27VoNGTJEUVFRSk5O1nvvvdfm8lu2bFFycrKioqJ02mmn6Yknnmi2zCuvvKKkpCRFRkYqKSlJGzdu7KrpAwCA41BIXwLbsGGDlixZorVr12rChAl68sknNX36dH366adKSEhotnxxcbF+8IMf6MYbb9T//M//aOvWrVq0aJH69eunK664QpK0bds2paWl6a677tLll1+ujRs36qqrrtL777+v888/v7tbBCBeTgNw7AnpFaCHHnpI8+fP14IFCzRy5EhlZGQoPj5e69ata3H5J554QgkJCcrIyNDIkSO1YMEC3XDDDfrd735nL5ORkaGpU6dq2bJlGjFihJYtW6ZLLrlEGRkZ3dQVAAA41oXsClB9fb3y8vJ02223BYynpqYqNze3xXW2bdum1NTUgLFp06bp6aefVkNDg8LDw7Vt2zYtXbq02TIEIODEwNUkAMEQsgBUUVEhn8+n2NjYgPHY2FiVl5e3uE55eXmLy3u9XlVUVGjAgAGtLtPaNiWprq5OdXV19vdVVVWSpH//+9+SJJ/PJ0lyu90BtdfrlWEYdu1yueRyuVqtGxoa5Ha7ZdZ9owi3pQafZMloUksRbqneJxmSwt1Svc+QIcuuDxw4IJ/Pp/DwcJmmKXfDITWYhlyGJbchNZiG3IYl1xHqMMPSv//974CezLpvFO6yZFqSzzKa1T5LMhtrUzJlKMJlyWtK1dXVamhoUFhYmAzDUJj3ULt7cslSmEuqN/11ZWWlwsLCZJqm/1H3Tbt7kiRvk7lXV1cHHCd3w6F292QGHBtD+/fvV1iY/1cmzHuoQz25XQo4TtXV1e3qr6WeDj82lZWV9jnmbjjUoZ4OP/caz32v1yuz7psO9dR07tXV1fL5fLIsy39O1B/qUE+t9ef1eqW6Qx3qqelxqq6ulmVZ8nq9Cg8P19kr3253T4fXf/vtlIB/F5Lvym53T4efe3m/SQ34N6Kxbvx9aqwbj03TOjw8PKAn0zQD/o0wTTPg9+nwuulxau3fumD8u0dP9NQdPVVWVkqSLMtq4b/2gUL+Z/CGYQR8b1lWs7EjLX/4eEe3uXr1at15553NxgcPHtzqOqHSKyN42+r9cPC2FZMRvG31CeK2YoLYY9+M4G0r5qHgbSuY++uUB4OznWO1v5gg9SdJvTOCt61eQdxfAKSamhrFxMS0uUzIAlDfvn3ldrubXZnZt29fsys4jeLi4lpcPiwsTH369Glzmda2KUnLli1Tenq6/b1pmqqsrFSfPn3aDE6S//8o4+PjVVpaqujo6DaXRXCx70OHfR867PvQYL+HTkf2vWVZqqmp0cCBA4+43ZAFoIiICCUnJys7O1uXX365PZ6dna1Zs2a1uE5KSor+7//+L2AsKytL48aNU3h4uL1MdnZ2wH1AWVlZGj9+fKtziYyMVGRkZMDYKaec0qF+oqOj+aUIEfZ96LDvQ4d9Hxrs99Bp774/0pWfRiF9CSw9PV1z5szRuHHjlJKSoqeeekolJSVauHChJP+Vma+++krPPfecJGnhwoV67LHHlJ6erhtvvFHbtm3T008/rfXr19vbvPXWWzVp0iTdf//9mjVrll577TW98847ev/990PSIwAAOPaENAClpaVp//79WrVqlcrKyjRq1Ch5PB4lJiZKksrKylRSUmIvP2TIEHk8Hi1dulSPP/64Bg4cqDVr1tjvASRJ48eP10svvaQVK1boN7/5jYYOHaoNGzbwHkAAAMAW8pugFy1apEWLFrX4XGZmZrOxyZMna+fOnW1u88orr9SVV14ZjOkdUWRkpFauXNnsJTR0PfZ96LDvQ4d9Hxrs99Dpqn1vWO35WzEAAIATSMg/CwwAAKC7EYAAAIDjEIAAAIDjEIAAAIDjEICO0tq1azVkyBBFRUUpOTlZ7733XqindMK74447ZBhGwCMuLi7U0zoh5eTk6LLLLtPAgQNlGIb+9Kc/BTxvWZbuuOMODRw4UD169NCFF16oTz75JDSTPYEcab9fd911zX4HLrjggtBM9gSyevVqnXvuuerZs6f69++v2bNn6/PPPw9YhnO+a7Rn3wf7vCcAHYUNGzZoyZIlWr58ufLz8zVx4kRNnz494L2L0DXOPPNMlZWV2Y+PPvoo1FM6IR06dEjnnHOOHnvssRaff+CBB/TQQw/pscce04cffqi4uDhNnTpVNTU13TzTE8uR9rskXXrppQG/Ax6PpxtneGLasmWLFi9erA8++EDZ2dnyer1KTU3VoUOH7GU457tGe/a9FOTz3kKnnXfeedbChQsDxkaMGGHddtttIZqRM6xcudI655xzQj0Nx5Fkbdy40f7eNE0rLi7Ouu++++yx2tpaKyYmxnriiSdCMMMT0+H73bIsa968edasWbNCMh8n2bdvnyXJ2rJli2VZnPPd6fB9b1nBP++5AtRJ9fX1ysvLU2pqasB4amqqcnNzQzQr5ygqKtLAgQM1ZMgQXX311dq9e3eop+Q4xcXFKi8vD/gdiIyM1OTJk/kd6Abvvvuu+vfvr2HDhunGG2/Uvn37Qj2lE05VVZUkqXfv3pI457vT4fu+UTDPewJQJ1VUVMjn8zX7lPnY2Nhmn0aP4Dr//PP13HPPadOmTfr973+v8vJyjR8/Xvv37w/11Byl8Tznd6D7TZ8+XS+88II2b96sBx98UB9++KEuvvhi1dXVhXpqJwzLspSenq7vf//7GjVqlCTO+e7S0r6Xgn/eh/yjMI53hmEEfG9ZVrMxBNf06dPt+qyzzlJKSoqGDh2qP/7xj0pPTw/hzJyJ34Hul5aWZtejRo3SuHHjlJiYqDfffFM/+tGPQjizE8fNN9+sXbt2tfhB2pzzXau1fR/s854rQJ3Ut29fud3uZql/3759zf7vAF3re9/7ns466ywVFRWFeiqO0viXd/wOhN6AAQOUmJjI70CQ3HLLLXr99df1l7/8Raeeeqo9zjnf9Vrb9y052vOeANRJERERSk5OVnZ2dsB4dna2xo8fH6JZOVNdXZ0KCws1YMCAUE/FUYYMGaK4uLiA34H6+npt2bKF34Futn//fpWWlvI7cJQsy9LNN9+sV199VZs3b9aQIUMCnuec7zpH2vctOdrznpfAjkJ6errmzJmjcePGKSUlRU899ZRKSkq0cOHCUE/thPbLX/5Sl112mRISErRv3z7dfffdqq6u1rx580I9tRPOwYMH9Y9//MP+vri4WAUFBerdu7cSEhK0ZMkS3XvvvTrjjDN0xhln6N5779VJJ52ka665JoSzPv61td979+6tO+64Q1dccYUGDBigPXv26Pbbb1ffvn11+eWXh3DWx7/FixfrxRdf1GuvvaaePXvaV3piYmLUo0cPGYbBOd9FjrTvDx48GPzzPmh/T+ZQjz/+uJWYmGhFRERYY8eODfiTPXSNtLQ0a8CAAVZ4eLg1cOBA60c/+pH1ySefhHpaJ6S//OUvlqRmj3nz5lmW5f+z4JUrV1pxcXFWZGSkNWnSJOujjz4K7aRPAG3t92+++cZKTU21+vXrZ4WHh1sJCQnWvHnzrJKSklBP+7jX0j6XZD377LP2MpzzXeNI+74rznvjPz8YAADAMbgHCAAAOA4BCAAAOA4BCAAAOA4BCAAAOA4BCAAAOA4BCAAAOA4BCAAAOA4BCAAAOA4BCEC3Ky0t1fz58zVw4EBFREQoMTFRt956q/bv39/ubezZs0eGYaigoKDrJgrghEUAAtCtdu/erXHjxunvf/+71q9fr3/84x964okn9Oc//1kpKSmqrKwM9RQBOAABCEC3Wrx4sSIiIpSVlaXJkycrISFB06dP1zvvvKOvvvpKy5cvlyQZhqE//elPAeuecsopyszMlCT706LHjBkjwzB04YUX2ss988wzOvPMMxUZGakBAwbo5ptvtp8rKSnRrFmzdPLJJys6OlpXXXWVvv76a/v5O+64Q6NHj9YzzzyjhIQEnXzyyfrZz34mn8+nBx54QHFxcerfv7/uueeegLlVVVXppptuUv/+/RUdHa2LL75Yf/vb3+zn//a3v+miiy5Sz549FR0dreTkZO3YsSMYuxRAJxCAAHSbyspKbdq0SYsWLVKPHj0CnouLi9O1116rDRs2qD0fUbh9+3ZJ0jvvvKOysjK9+uqrkqR169Zp8eLFuummm/TRRx/p9ddf1+mnny5JsixLs2fPVmVlpbZs2aLs7Gz985//VFpaWsC2//nPf+qtt97S22+/rfXr1+uZZ57RjBkz9OWXX2rLli26//77tWLFCn3wwQf2dmfMmKHy8nJ5PB7l5eVp7NixuuSSS+wrWtdee61OPfVUffjhh8rLy9Ntt92m8PDwo9uhADrv6D/DFQDa54MPPrAkWRs3bmzx+YceesiSZH399dctLhcTE2N/OnRxcbElycrPzw9YZuDAgdby5ctb3H5WVpbldrsDPkH6k08+sSRZ27dvtyzLslauXGmddNJJVnV1tb3MtGnTrMGDB1s+n88eGz58uLV69WrLsizrz3/+sxUdHW3V1tYG/LyhQ4daTz75pGVZltWzZ08rMzOz5R0DoNtxBQjAMcP6z5UfwzA6tf6+ffu0d+9eXXLJJS0+X1hYqPj4eMXHx9tjSUlJOuWUU1RYWGiPDR48WD179rS/j42NVVJSklwuV8DYvn37JEl5eXk6ePCg+vTpo5NPPtl+FBcX65///KckKT09XQsWLNCUKVN033332eMAQoMABKDbnH766TIMQ59++mmLz3/22Wfq1auX+vbtK8Mwmr0U1tDQ0Ob2D39Z7XCWZbUYrg4fP/ylKcMwWhwzTVOSZJqmBgwYoIKCgoDH559/rl/96leS/PcWffLJJ5oxY4Y2b96spKQkbdy4sc35Aug6BCAA3aZPnz6aOnWq1q5dq2+//TbgufLycr3wwgtKS0uTYRjq16+fysrK7OeLior0zTff2N9HRERIknw+nz3Ws2dPDR48WH/+859b/PlJSUkqKSlRaWmpPfbpp5+qqqpKI0eO7HRfY8eOVXl5ucLCwnT66acHPPr27WsvN2zYMC1dulRZWVn60Y9+pGeffbbTPxPA0SEAAehWjz32mOrq6jRt2jTl5OSotLRUb7/9tqZOnapBgwbZf1118cUX67HHHtPOnTu1Y8cOLVy4MOAqTP/+/dWjRw+9/fbb+vrrr1VVVSXJf6XlwQcf1Jo1a1RUVKSdO3fq0UcflSRNmTJFZ599tq699lrt3LlT27dv19y5czV58mSNGzeu0z1NmTJFKSkpmj17tjZt2qQ9e/YoNzdXK1as0I4dO/Ttt9/q5ptv1rvvvqsvvvhCW7du1YcffnhUoQvA0SEAAehWZ5xxhnbs2KGhQ4cqLS1NQ4cO1U033aSLLrpI27ZtU+/evSVJDz74oOLj4zVp0iRdc801+uUvf6mTTjrJ3k5YWJjWrFmjJ598UgMHDtSsWbMkSfPmzVNGRobWrl2rM888UzNnzlRRUZGk7/60vlevXpo0aZKmTJmi0047TRs2bDiqngzDkMfj0aRJk3TDDTdo2LBhuvrqq7Vnzx7FxsbK7XZr//79mjt3roYNG6arrrpK06dP15133nlUPxdA5xnW4S+yAwAAnOC4AgQAAByHAAQAAByHAAQAAByHAAQAAByHAAQAAByHAAQAAByHAAQAAByHAAQAAByHAAQAAByHAAQAAByHAAQAAByHAAQAABzn/wPmGGx6h1bijAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "P6 = monte_carlo(roll_1, 6)\n",
    "P_barchart(P6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that 1/6th of the time the game ends on the first roll, but 2 or 3 rolls are more common, and games of 20 or more rolls do occur. We can compute the average, or expected number of rolls (by convention `E` is used for the expected value):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.287199"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def E(P: ProbDist) -> float:\n",
    "    \"\"\"The expected value of a probability distribution.\"\"\"\n",
    "    return sum(P[x] * x for x in P)\n",
    "\n",
    "E(P6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An average game takes about 3.28 rolls. \n",
    "\n",
    "What is the probability that the first player wins? That's the sum of the probabilities of the outcomes that have an odd number of rolls:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.523025"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def first_player_wins(P: ProbDist) -> Probability:\n",
    "    \"\"\"The probability that the first player wins.\"\"\"\n",
    "    return sum(P[x] for x in P if x % 2 == 1)\n",
    "\n",
    "first_player_wins(P6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A simulation with a million rolls should give you about 3 digits of accuracy.  If that degree of (in)accuracy is what you want, then we're done. \n",
    "\n",
    "# Other Approaches to Roll-1\n",
    "\n",
    "The simulation was easy, but could we use a formula or a complete enumeration?\n",
    "\n",
    "## Formula for Roll-1\n",
    "\n",
    "I can give a formula for  the special  case of a 2-sided die. This is essentially a series of coin flips: on each turn there is a ½ chance of a 1 and ½ chance of a 2. The first player wins if there is an even number of 2s (zero or more) followed by a 1.  So the total probability for the first player to win is the following infinite [geometric series](https://en.wikipedia.org/wiki/Geometric_series), for which there is a [formula](https://en.wikipedia.org/wiki/Geometric_series#Formula):\n",
    "\n",
    "\n",
    "$$½ + ½^3 + ½^5 + ½^7 + \\ldots = a / (1 - r) = ½ / (1 - ½^2) = 2/3.$$\n",
    "\n",
    "But if we start with a die with more than 2 sides then I don't know of a concise mathematical formula. \n",
    "\n",
    "## Enumerating Probabilities for Roll-1\n",
    "\n",
    "If we could enumerate every possible game of Roll-1, we could calculate the *exact* win percentage and expected number of rolls. The bad news is that the game is **unbounded**: any number of rolls is possible, so no finite enumeration can capture every possible game. (This is sometimes called a **non-well-founded game**.)  The good news is that games with a huge number of rolls are *very unlikely*. Thus, we can enumerate every possibility for a finite, truncated version of the game, and that should give results that are very close to the full game. We will consider all possible games up to a depth of, say, 100 rolls. For a roll there are *d* possibilities. If a player rolls a 1, they win, with certainty (that's the 1.0 in the `return` statement below). For each of the die rolls from 2 to *d*, the probability that the first player wins is 1 minus the probability that the other player wins (when they become the first player in the remainder of the game). The final probability is the mean of these *d* outcomes. If we reach *depth* rolls, we truncate the game and assume that both players are equally likely to win. (This is not quite right, but since a game with over 100 rolls happens vary rarely, the overall error will be small.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "from functools import lru_cache\n",
    "\n",
    "@lru_cache(None)\n",
    "def roll_1_first_player_wins(d, depth=100) -> Probability:\n",
    "    \"\"\"The probability that the first player wins a game of Roll-1.\"\"\"\n",
    "    if depth == 0:\n",
    "        return 1/2 # A guess\n",
    "    else: # First player wins by rolling a 1, or wins with probability (1.0 - the_other_player_wins) for other rolls.\n",
    "        second_player_loses = [1.0 - roll_1_first_player_wins(d1, depth - 1) for d1 in range(2, d + 1)]\n",
    "        return mean([1.0, *second_player_loses])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are the win probabilitues for 2- and 6-sided dice:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6666666666666666"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "roll_1_first_player_wins(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5238095238095237"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "roll_1_first_player_wins(6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What about computing the average length of a game? The function `roll_1_average_length` does this, using the same strategy as above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "@lru_cache(None)\n",
    "def roll_1_average_length(d, depth=100) -> float:\n",
    "    \"\"\"The average length of a game of Roll-1 (truncated to `depth` moves).\"\"\"\n",
    "    if depth == 0:\n",
    "        return 1 # A guess\n",
    "    else:\n",
    "        lengths = [1 + roll_1_average_length(d1, depth - 1) for d1 in range(2, d + 1)]\n",
    "        return mean([1, *lengths])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.0"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "roll_1_average_length(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.283333333333333"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "roll_1_average_length(6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try some other values of *d*, and compare the simulation approach with the exact (albeit truncated) approach:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average game length and win probability for various-sided dice, with simulated/exact calculations:\n",
      "\n",
      "| sides | sim/exact length | sim/exact win % |\n",
      "|     2 | 1.9996 / 2.0000  | 66.62% / 66.67% |\n",
      "|     3 | 2.4994 / 2.5000  | 58.38% / 58.33% |\n",
      "|     4 | 2.8334 / 2.8333  | 54.97% / 55.00% |\n",
      "|     5 | 3.0863 / 3.0833  | 53.27% / 53.33% |\n",
      "|     6 | 3.2836 / 3.2833  | 52.33% / 52.38% |\n",
      "|     8 | 3.5968 / 3.5929  | 51.39% / 51.39% |\n",
      "|    10 | 3.8279 / 3.8290  | 50.90% / 50.91% |\n",
      "|    12 | 4.0177 / 4.0199  | 50.62% / 50.64% |\n",
      "|    20 | 4.5431 / 4.5477  | 50.23% / 50.24% |\n",
      "|    32 | 5.0285 / 5.0272  | 50.13% / 50.09% |\n",
      "|    60 | 5.6703 / 5.6632  | 50.09% / 50.03% |\n",
      "|   100 | 6.1799 / 6.1774  | 50.02% / 50.01% |\n"
     ]
    }
   ],
   "source": [
    "def show(sides=[2, 3, 4, 5, 6, 8, 10, 12, 20, 32, 60, 100]) -> None:\n",
    "    \"\"\"Print information on simulated/exact game lengths and win percentages for Roll-1.\"\"\"\n",
    "    print('Average game length and win probability for various-sided dice, with simulated/exact calculations:\\n')\n",
    "    print('| sides | sim/exact length | sim/exact win % |')\n",
    "    for d in sides:\n",
    "        P    = monte_carlo(roll_1, d)\n",
    "        avg1 = E(P)\n",
    "        win1 = first_player_wins(P)\n",
    "        avg2 = roll_1_average_length(d)\n",
    "        win2 = roll_1_first_player_wins(d)\n",
    "        print(f'| {d:5} | {avg1:6.4f} / {avg2:6.4f}  | {win1:5.2%} / {win2:5.2%} |')\n",
    "\n",
    "show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "The simulation and exact calculations yield similar results, usually to 3 decimal places. But how exact is the \"exact\" calculation? There are two potential sources of error:\n",
    "- Python uses 64-bit floating point numbers, with 52-bit mantissas (the part after the decimal point). That  allows for no more than 16 digits of accuracy.\n",
    "- We truncate the search after *depth* rolls, taking a guess at the remaining computation.\n",
    "\n",
    "We can measure the error lost by truncation by comparing the difference between results at  *depth* and results at deeper depths."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def diff(fn: Callable, d: int, depth=100, deeper=120) -> List[float]:\n",
    "    \"\"\"The maximum difference between fn(d, depth) and fn(d, depth2) for depth2 up to `deeper`.\"\"\"\n",
    "    return max(abs(fn(d, depth) - fn(d, depth2)) for depth2 in range(depth, deeper + 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1.1102230246251565e-16,\n",
       " 1.1102230246251565e-16,\n",
       " 1.1102230246251565e-16,\n",
       " 0.0,\n",
       " 1.1102230246251565e-16,\n",
       " 1.1102230246251565e-16,\n",
       " 1.1102230246251565e-16,\n",
       " 1.1102230246251565e-16,\n",
       " 1.1102230246251565e-16,\n",
       " 0.0,\n",
       " 1.1102230246251565e-16,\n",
       " 1.1102230246251565e-16,\n",
       " 0.0,\n",
       " 0.0,\n",
       " 0.0,\n",
       " 0.0,\n",
       " 1.1102230246251565e-16,\n",
       " 0.0,\n",
       " 0.0,\n",
       " 0.0,\n",
       " 0.0,\n",
       " 0.0]"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[diff(roll_1_first_player_wins, d) for d in range(2, 24)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.0,\n",
       " 0.0,\n",
       " 0.0,\n",
       " 0.0,\n",
       " 0.0,\n",
       " 0.0,\n",
       " 0.0,\n",
       " 0.0,\n",
       " 0.0,\n",
       " 0.0,\n",
       " 0.0,\n",
       " 0.0,\n",
       " 0.0,\n",
       " 0.0,\n",
       " 0.0,\n",
       " 0.0,\n",
       " 0.0,\n",
       " 0.0,\n",
       " 0.0,\n",
       " 0.0,\n",
       " 0.0,\n",
       " 0.0]"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[diff(roll_1_average_length, d) for d in range(2, 24)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The number 1.1102230246251565e-16 is 2<sup>-53</sup>, the smallest possible difference between two 64-bit floating point numbers (also known as the \"unit in the last place\" or ULP). So truncating at depth 100 leads to little or no error."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Continuous Sample Spaces\n",
    "\n",
    "Everything up to here has been about discrete sample spaces, where we can enumerate all the possible outcomes.\n",
    "\n",
    "But what about *continuous* sample spaces, such as the space of real numbers? The principles are the same: probability is still the ratio of the favorable cases to all the cases. It is still easy to do a simulation. But for an exact computation, instead of *enumerating* all the discrete cases, we have to (in general) compute integrals. \n",
    "Here we will cover a simple continuous example, which we first solve approximately by simulation, and then exactly by calculation.\n",
    "\n",
    "## The Hot New Game Show Problem: Simulation\n",
    "\n",
    "[Oliver Roeder](https://fivethirtyeight.com/contributors/oliver-roeder/) posed [this problem](http://fivethirtyeight.com/features/can-you-win-this-hot-new-game-show/) in the 538 *Riddler* blog:\n",
    "\n",
    ">Two players go on a hot new game show called *Higher Number Wins.* The two go into separate booths, and each presses a button, and a random number between zero and one appears on a screen. (At this point, neither knows the other’s number, but they do know the numbers are chosen from a standard uniform distribution.) They can choose to keep that first number, or to press the button again to discard the first number and get a second random number, which they must keep. Then, they come out of their booths and see the final number for each player on the wall. The lavish grand prize — a case full of gold bullion — is awarded to the player who kept the higher number. Which number is the optimal cutoff for players to discard their first number and choose another? Put another way, within which range should they choose to keep the first number, and within which range should they reject it and try their luck with a second number?\n",
    "\n",
    "We'll use this notation:\n",
    "- **A**, **B**: the two players.\n",
    "- *A*, *B*: the cutoff values they choose: the lower bound of the range of first numbers they will accept.\n",
    "- *a*, *b*: the actual random numbers that appear on the screen.\n",
    "\n",
    "For example, if player **A** chooses a cutoff of *A* = 0.6, that means that **A** would accept any first number greater than 0.6, and reject any number below that cutoff. The question is: What cutoff, *A*, should player **A** choose to maximize the chance of winning, that is, maximize P(*a* > *b*)? (Since the situation is symmetric, the cutoff value *B* will be the same.)\n",
    "\n",
    "First, simulate a single player with a given cutoff playing the game once and getting a number:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "number = random.random\n",
    "\n",
    "def get_number(cutoff) -> float:\n",
    "    \"\"\"Simulate a single player with a given cutoff receiving their number.\"\"\"\n",
    "    first = number()\n",
    "    return first if first > cutoff else number()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compare two strategies: averaged over many simulated trials, how often does a player with cutoff *A* defeat a player with cutoff *B*?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Pwin_sim(A, B, trials=100_000) -> Probability:\n",
    "    \"The probability that cutoff A wins against cutoff B.\"\n",
    "    return mean(get_number(A) > get_number(B) for _ in range(trials))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.52916"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Pwin_sim(0.6, 0.8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A cutoff of 0.6 appears to be better than a cutoff of 0.8, winning about 53% of the time.\n",
    "\n",
    "Now define a function, `top`, that considers a list of possible cutoffs, estimates the win probability for each cutoff playing against each other cutoff, and returns a list with the `N` top cutoffs (the ones that defeated the most number of opponent cutoffs), and the number of opponents they defeat: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "def arange(start: float, stop: float, step: float) -> List[float]:\n",
    "    \"\"\"Like numpy.arange, a list of numbers from start to stop by step.\"\"\"\n",
    "    length = int((stop - start) / step)\n",
    "    return [start + i * step for i in range(length)]\n",
    "\n",
    "def top(N: int, Pwin_fn: Callable, cutoffs: List[float]) -> List[Tuple[float, int]]:\n",
    "    \"Return the N best cutoffs and the number of opponent cutoffs they beat.\"\n",
    "    winners = Counter(A if Pwin_fn(A, B) > 0.5 else B\n",
    "                      for (A, B) in combinations(cutoffs, 2))\n",
    "    return winners.most_common(N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(0.58, 46),\n",
       " (0.59, 46),\n",
       " (0.6, 46),\n",
       " (0.63, 46),\n",
       " (0.62, 45),\n",
       " (0.61, 43),\n",
       " (0.56, 42),\n",
       " (0.5700000000000001, 42),\n",
       " (0.55, 41),\n",
       " (0.64, 41)]"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "top(10, Pwin_sim, arange(0.50, 1.00, .01))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This suggests that the best cutoff is somewhere around 0.60, but the results do not seem definitive. We could run more simulations, but that would take a lot of computation time.\n",
    "\n",
    "## The Hot New Game Show Problem: Exact Calculation\n",
    "\n",
    "More promising is the possibility of replacing the simulation with an exact calculation. But first let's solve a simpler problem: assume that both players **A** and **B** have chosen a cutoff, and have each received a first number that is above the cutoff. What is the probability that **A** gets the higher number? We'll call this `Phigher(A, B)`. We can think of this as a continuous two-dimensional sample space of points in the (*a*, *b*) plane, where *a* ranges from the cutoff *A* to 1 and *b* ranges from the cutoff *B* to 1. Here is a diagram of that two-dimensional sample space, with the cutoffs *A*=0.5 and *B*=0.6:\n",
    "\n",
    "<img src=\"http://norvig.com/ipython/probability2da.jpg\" width=413>\n",
    "\n",
    "The total area of the sample space is 0.5 &times; 0.4 = 0.20, and in general it is (1 - *A*) &times; (1 - *B*). What about the  cases where **A** beats **B**? That corresponds to the shaded triangle below:\n",
    "\n",
    "<img src=\"http://norvig.com/ipython/probability2d.jpg\" width=413>\n",
    "\n",
    "The area of a triangle is 1/2 the base times the height, or in this case, 0.4<sup>2</sup> / 2 = 0.08, and in general, (1 - *B*)<sup>2</sup> / 2. So in general we have:\n",
    "\n",
    "    Phigher(A, B) = favorable / total\n",
    "    favorable     = ((1 - B) ** 2) / 2 \n",
    "    total         = (1 - A) * (1 - B)\n",
    "    Phigher(A, B) = (((1 - B) ** 2) / 2) / ((1 - A) * (1 - B))\n",
    "    Phigher(A, B) = (1 - B) / (2 * (1 - A))\n",
    "    \n",
    "And in this specific case we have:\n",
    "\n",
    "     A = 0.5; B = 0.6\n",
    "     favorable         = 0.4 ** 2 / 2 = 0.08\n",
    "     total             = 0.5 * 0.4    = 0.20\n",
    "     Phigher(0.5, 0.6) = 0.08 / 0.20  = 0.4\n",
    "\n",
    "But note that this only works when the cutoff *A* &le; *B*; when *A* > *B*, we need to reverse things. That gives us the code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Phigher(A, B):\n",
    "    \"\"\"Probability that a sample from [A..1] is higher than one from [B..1].\"\"\"\n",
    "    if A <= B:\n",
    "        return (1 - B) / (2 * (1 - A))\n",
    "    else:\n",
    "        return 1 - Phigher(B, A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert Phigher(0.5, 0.6) == 0.4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're now ready to tackle the full game. There are four cases to consider, depending on whether **A** and **B** gets a first number that is above or below their cutoff choices:\n",
    "\n",
    "| first *a* | first *b* | P(*a*, *b*) | P(A wins &vert; *a*, *b*) |  Comment |\n",
    "|:-----:|:-----:| ----------- | ------------- | ------------ |\n",
    "| *a* > *A* | *b* > *B* | (1 - *A*) &middot; (1 - *B*) | Phigher(*A*, *B*) | Both above cutoff; both keep first numbers |\n",
    "| *a* < *A* | *b* < *B* | *A* &middot; *B* | Phigher(0, 0)  | Both below cutoff, both get new numbers from [0..1] |\n",
    "| *a* > *A* | *b* < *B* | (1 - *A*) &middot;  *B* | Phigher(*A*, 0) | **A** keeps number; **B** gets new number from [0..1] |\n",
    "| *a* < *A* | *b* > *B* | *A* &middot; (1 - *B*) | Phigher(0, *B*) | **A** gets new number from [0..1]; **B** keeps number |\n",
    "\n",
    "For example, the first row of this table says that the event of both first numbers being above their respective cutoffs has probability (1 - *A*) &middot; (1 - *B*), and if this does occur, then the probability of **A** winning is Phigher(*A*, *B*).\n",
    "We're ready to replace the old simulation-based `Pwin_sim` with a new exact version, `Pwin_exact`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Pwin_exact(A, B):\n",
    "    \"With what probability does cutoff A win against cutoff B?\"\n",
    "    return ((1-A) * (1-B) * Phigher(A, B)  # both above cutoff\n",
    "            + A * B       * Phigher(0, 0)  # both below cutoff\n",
    "            + (1-A) * B   * Phigher(A, 0)  # A above, B below\n",
    "            + A * (1-B)   * Phigher(0, B)) # A below, B above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.528"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Pwin_exact(0.6, 0.8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the 0.6 cutoff beats 0.8 52.8% of the time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Pwin_exact` relies on a lot of algebra. Let's define a few tests to check for obvious errors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'ok'"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def test():\n",
    "    assert Phigher(0.5, 0.5) == Phigher(0.75, 0.75) == Phigher(0, 0) == 0.5\n",
    "    assert Pwin_exact(0.5, 0.5) == Pwin_exact(0.75, 0.75) == 0.5\n",
    "    assert Phigher(.6, .5) == 0.6\n",
    "    assert Phigher(.5, .6) == 0.4\n",
    "    return 'ok'\n",
    "\n",
    "test()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's repeat the `top` calculation with `Pwin_exact`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(0.62, 49),\n",
       " (0.61, 48),\n",
       " (0.6, 47),\n",
       " (0.59, 46),\n",
       " (0.63, 45),\n",
       " (0.58, 44),\n",
       " (0.5700000000000001, 43),\n",
       " (0.64, 42),\n",
       " (0.56, 41),\n",
       " (0.55, 40)]"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "top(10, Pwin_exact, arange(0.5, 1.0, 0.01))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see here that of the 50 candidate cutoffs, 0.62 defeats all the other 49, and 0.61 defeats 48 (all except 0.62). The great thing about the exact calculation code is that it runs fast, regardless of how much accuracy we want. We can zero in on those top ten cutoffs, 0.55 to 0.64, and go to 4 digits:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(0.6180000000000001, 898),\n",
       " (0.6179, 897),\n",
       " (0.6181000000000001, 896),\n",
       " (0.6178, 895),\n",
       " (0.6177, 894),\n",
       " (0.6176, 893),\n",
       " (0.6182000000000001, 892),\n",
       " (0.6175, 891),\n",
       " (0.6174000000000001, 890),\n",
       " (0.6183000000000001, 889)]"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "top(10, Pwin_exact, arange(0.55, 0.64, 0.0001))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This says 0.6180 is best. Let's go to 6 digits:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(0.618034, 899),\n",
       " (0.6180329999999999, 898),\n",
       " (0.6180319999999999, 897),\n",
       " (0.618035, 896),\n",
       " (0.618031, 895),\n",
       " (0.61803, 894),\n",
       " (0.6180289999999999, 893),\n",
       " (0.6180359999999999, 892),\n",
       " (0.6180279999999999, 891),\n",
       " (0.618027, 890)]"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "top(10, Pwin_exact, arange(0.6174, 0.6183, 0.000001))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now 0.618034 is best. I'll leave it up to you to explore whether the exact answer is actually a [well-known irrational number](https://en.wikipedia.org/wiki/Golden_ratio).\n",
    "\n",
    "So far we've been hunting around in the space of possible cutoffs, trying to find the best one, but we can make a stronger argument. Suppose the game is changed so that player **A** has to declare a cutoff first, and then player **B** gets to respond with their cutoff, with full knowledge of **A**'s choice. In other words, what cutoff should **A** choose to **maximize** `Pwin_exact(A, B)`, given that **B** is going to take that knowledge and pick a cutoff that **minimizes** `Pwin_exact(A, B)`? \n",
    "\n",
    "We can consider a wide range of cutoffs, with some that zero in on the interval around 0.618034."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "cutoffs = list(arange(0.000, 1.000, 0.001)) + list(arange(0.618030, 0.618040, 0.000001))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I'll use this helper function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Pwin(A, B): return (Pwin_exact(A, B), f'probability of winning for A with A={A}, B={B}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is the final winning percentage for **A** along with the best choice of cutoffs *A* and *B*, when **A** declares their cutoff first:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.5, 'probability of winning for A with A=0.618034, B=0.618034')"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "max(min(Pwin(A, B) for B in cutoffs)\n",
    "    for A in cutoffs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here it is the other way around, when **B** has to declare their cutoff first:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.5, 'probability of winning for A with A=0.618034, B=0.618034')"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "min(max(Pwin(A, B) for A in cutoffs)\n",
    "    for B in cutoffs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In both cases, the rational choice for both players in a cutoff of 0.618034. This is a *stable equilibrium*: consider fixing *B* = 0.618034, and notice that if *A* changes to any other cutoff value, player **A** does worse. Similarly, if we fix *A* = 0.618034, then if *B* changes to another value, player **B** does worse. So neither player will want to move.\n",
    "\n",
    "The moral for continuous spaces is the same as for discrete spaces: be careful about defining your sample space; measure carefully, and let your code take care of the rest."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
